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Notations and symbols 



e Element of 

C Subset of 

= IdentiacUy equals 

~ Essentially equal to or equivalent 

f Hermitian conjugate 

/+, ^4+ Generalised solution and generalised inverse respectively 

/' First derivative of / 

(/) Mean value of / 

X — {x^) Vector in Euclidean n-dimensional space R". The com- 
ponents are labelled by latin letters (j = 1,2,3) 

X = -. — j- Unit vector in the direction of x 
\x\ 

X = (x'*) 4- Vector in 4-dimensional Minkowski space. The com- 
ponents are labelled by greek letters (// = 0, 1, 2, 3) 

(T'x, dx Volume element in Euclidean n-dimensional space 

d'^x — dx^dx^dx'^dx^ Volume element in 4-dimensional Minkowski space 

/ Unit operator or matrix 

n Closure of the domain VL 

do, Boundary of the domain O 

0{E) Terms of order E 

N{A) Null-space of operator A 

Ti-^A) Range of operator A 

T>[A) Domain of operator A 

A* Adjoint of A 

A'^ Transposed of A 

TrA Trace of A 

I . I Absolute value 

Norm defined on the space X 



Re X Real part of x 

Im X Imaginary part of x 

P Principal part 

Summation 

[0\ Mass dimesion of operator O, i.e., [0\ = M'^C^) 

T[f]{x) Fourier transform of / 

(/> 9)x Scalar product defined on the space X 

A Laplace operator 

V Nabla operator 

5{x) Dirac (5-function {—oo < x < oo); 6{x) = for a; 7^ 0, 

J^^dxS{x) = 1 

9{x) Heaviside-function (—00 < a; < 00); 9{x) = for a; < 0, 

e{x) = 1 for X > 

g^, Quark fields. Flavours are labelled by latin letters, here 

i, and colours by greek letters, here a 

C^" Gluon field strength tensor 

7^^,75 Dirac matrices 

A" Gell-Mann colour matrices 

fabc Structure constant of SU{3) 

: A ■ B ■ ... : Normal ordering of the operators A,B, ... 



Abbreviations 



MC 


Monte Carlo 


LO 


Leading order 


NLO 


Next-to-leading order 


SVD 


Singular value decomposition 


QCD 


Quantum chromo dynamics 


QED 


Quantum electrodynamics 


OPE 


Operator product expansion 


CKM 


Cabibbo-Kobayashi-Maskawa 


PCAC 


Partial conservation of axial-vector current 


EIT 


Electrical impedance tomography 


FEM 


Finite element method 


CL 


Confidence level 


LS 


Least-squares 


CR 


Confidence region 


p.d.f. 


probability distribution function 
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Introduction 



We call two problems inverse to each other if the formulation of each of them 
requires full or partial knowledge of the other. By this definition, it is obviously 
arbitrary which of the two problems we call the direct and which we call the inverse 
problem. But usually, one of the problems has been studied earlier and, perhaps, in 
more detail. This one is usually called the direct problem, whereas the other is the 
inverse problem. However, there is often another, more important difference between 
these two problems. Hadamard (see |Ha23] ) introduced the concept of a well-posed 
problem, originating from the philosophy that the mathematical model of a physical 
problem has to have the properties of uniqueness, existence, and stability of the 
solution. If one of the properties fails to hold, he called the problem ill-posed. It turns 
out that many interesting and important inverse problems in science lead to ill-posed 
problems, while the corresponding direct problems are well-posed. Often, existence 
and uniqueness can be forced by enlarging or reducing the solution space (the space 
of "models"). For restoring stability, however, one has to change the topology of the 
space, which is in many cases impossible because of the presence of measurement 
errors. At first glance, it seems to be impossible to compute the solution of a problem 
numerically if the solution of the problem does not depend continuously on the 
data, i.e., for the case of ill-posed problems. Under additional a priori information 
about the solution, such as smoothness and bounds on the derivatives, however, it 
is possible to restore stability and construct efficient numerical algorithms. 

The thesis contains three main parts. The aim of the first part is to introduce 
the basic notations and difficulties encountered with ill-posed problems and then 
study the basic properties of regularisation methods for linear ill-posed problems. 
In the second and third part we aim to find stable solutions to two inverse problems 
arising in different fields of physics. 

The first is an inverse problem of quantum chromo dynamics (QCD). QCD is 
widely considered to be a good candidate for a theory of the strong interactions. 
Asymptotic freedom allows us to perform a perturbative treatment of strong inter- 
actions at short distances. Long distance behaviour is not fully understood: it is 
commonly believed that, due to the nontrivial structure of the physical vacuum, 
the perturbation expansion does not completely define the theory and that one has 
to add non-perturbative effects as well. In order to make a comparison with ex- 
periment possible, even in the resonance energy range, Shifman, Vainshtein and 
Zakharov |5VZ79a] have proposed to use the Operator Product Expansion (OPE) 
and to introduce the vacuum expectation values of the operators occurring in the 
OPE, the so called condensates, as phenomenological parameters. The knowledge of 



11 



12 



these parameters is useful for studying if one can indeed obtain a consistent descrip- 
tion of the low energy hadronic physics and get more insight into the properties of 
the QCD vacuum. It is therefore important not only to determine the values of the 
condensates from experimental data, but also the accuracy of the determination, i.e. 
to determine their allowed range. 

In QCD, perturbation expansions, together with some non-perturbative effects, 
allows us to approximate the two point functions of hadronic currents, i.e., vacuum 
expectation values of the time ordered product of two hadronic currents, in the 
distant space-like region in terms of a few parameters (the values of condensates). 
On the other hand, the discontinuity of these amplitudes in the time-like region 
is related to more directly measurable quantities. Although analyticity strongly 
correlates the values of the amplitudes in these two regions, the errors affecting 
both types of data make the correlation much looser. Any procedure aimed to build 
up acceptable amplitudes must take into account these errors in a reasonable way. 

There are several methods, generically called QCD sum rules, dealing with this 
problem of analytic extrapolation [BBSlaj IBLR85t ISVZ79a] . Most of them include 
the theoretical errors in the space-like region only at a qualitative level, and/or 
need (explicit and implicit) assumptions on the derivatives of the amplitudes. The 
application of fully controlled analytic extrapolation techniques should remedy these 
effects. There are a few methods of this sort, in which the error channels in the space- 
like region are defined through L^-norms |A5507| IC5S504] or L°°-norms |AM89 
lACMQOlirMQO] . 

The functional method we use [AS 5 07] allows us to extract within rather general 
assumptions the condensates from a comparison of the time-like experimental data 
with the asymptotic space-like results from theory. We will see that the price to be 
paid for the generality of assumptions is relatively large errors in the values of the 
extracted parameters. Although we do not claim that our method is superior to other 
approaches, we hope that our results lend additional confidence to the numerical 
results obtained with the help of methods based on QCD sum rules |BB81a[ IBBBlb 



[BiSTl rasl [BLR851 [D5881 [b05l [PCPSM ILNT841 [Ni95l [Ni95l [Ni96l [Ni98l [Ni02 
[Ni04l IRRY851 l5VZ79a| rYn99] . 

For the experimental data on the time-like region, one has more possibilities 
at hand. We have chosen to work with the final r data provided by the ALEPH 
collaboration [ ALEPH05] . because they have the smallest experimental errors. One 
may also ask why r-decays and not e^e~-annihilation? The answer lies in the fact 
that the physics of hadronic r-decays has been the subject of much progress in the 
last decade, both at the experimental and theoretical level. Somewhat unexpectedly, 
hadronic r-decays provide one of the most powerful testing grounds for QCD. This 
situation results from a number of favourable conditions: 

• The r lepton is heavy enough to decay into a variety of hadrons, with net 
strangeness and ±1; 

• r leptons are copiously produced in pairs at e+e^ colliders, leading to simple 
event topologies with little background; 
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• Experimental studies of r-decays could be performed with large data samples; 

• As a consequence, r-decays are experimentally known to great detail and their 
rates measured with high precision; 

• The theoretical description of hadronic r-decays is based on solid ground. 

The second inverse problem we aimed to solve is in the field of Electrical Im- 
pedance Tomography (EIT). EIT is a technology developed to image the electrical 
conductivity distribution of a conductive medium. It is of interest because of its low 
cost and also because the electrical conductivity gives direct information about the 
internal composition of the conductive medium. The technique works by performing 
simultaneous measurements of direct or alternating electric currents and voltages on 
the boundary of an object. These are the data used by an image reconstruction al- 
gorithm to determine the electrical conductivity distribution within the object. This 
problem is also called the inverse conductivity problem and has applications as a 
method of industrial, geophysical and medical imaging. 

It is well known that the inverse conductivity problem is a highly ill-posed, 
non-linear inverse problem and that the images produced are very sensitive to er- 
rors which can occur in practice. There has been much interest in determining the 
class of conductivity distributions that can be recovered from the boundary data, 
as well as in the development of related reconstruction algorithms. The interest 
in this problem has been generated by both difficult theoretical challenges and by 
the important medical, geophysical and industrial application of it. Much theore- 
tical work has been related to the approach of Calderon concerning the bijection 
between the the conductivity inside the region and the Neumann-to-Dirichlet oper- 
ator, which relates the distribution of the injected currents to the boundary values of 
the induced electrical potential [CaSOl [KV85t INa88l l5U88j . The reconstruction pro- 
cedures that have been proposed include a wide range of iterative methods based 
on formulating the inverse problem as a nonlinear optimisation problem. These 
techniques are quite demanding computationally, particularly when addressing the 
three-dimensional problem. This drawback has encouraged the search for recon- 
struction algorithms which reduce the computational demands either by using some 
a priori information e.g. [BHOOl IBHVOSt |CFMV98J or by developing non-iterative 
procedures. Some of these methods [AHQS07t iBROOl IBHVOSj use a factorisation 
approach while others are based on reformulating the inverse problem in terms of 
integral equations [CPSSOOl ICPSOSIITPSM] . 

One of the approaches presented here is based on reformulating the inverse prob- 
lem in terms of integral equations. A feature of this method is that many of the 
calculations involve analytical expressions containing the eigenfunctions of the ker- 
nel of these equations, the computational part being restricted to the introduction 
of the data, the numerical evaluation of some of the analytic formulae and the so- 
lution of a final integral equation. The method consists in the determination of 
Y{x) = —A(f){x), in the sense of a generalised solution of inverse problems, by a 
single measurement of the potential (pix), and its normal derivative on the bound- 
ary of the domain. The result Y{x) can either be used directly to obtain rough 
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information on the conductivity or may be processed further to determine the con- 
ductivity by solving a first order partial differential equation. This can be done in 
a straightforward way by using the method of characteristics. We have applied this 
method for a two-dimensional domain, the unit disc, with no a priori information. 
However, since the problem is ill-posed, one needs to regularise it, i.e., search for 
approximate solutions satisfying additional constraints suggested by the physics of 
the problem. In our case we have used as a regularisation algorithm the truncated 
singular value decomposition. Unfortunately, the information gained on Y{x), and 
the conductivity respectively, is restricted to their angular dependence (no radial 
information is present). One can hope that using some a priori information could 
improve the reconstruction. 

Also an algorithm based on linearisation will be discussed. The aim of the algo- 
rithm is to perform reconstructions based on real data measured by the tomograph 
constructed in collaboration with Dr. K.H. Georgi and N. Schuster. For this, a 
belt of electrodes was placed around the chest of a human volunteer, voltages ap- 
plied and currents measured. An important result of this algorithm is that the low 
conductivity (lungs) and high conductivity (heart) regions are well reconstructed. 
This is important for monitoring for lung problems such as accumulating fluid or a 
collapsed lung. 

The EIT research is still active. There are about 30 groups worldwide who 
are actively performing research and it is still seen as an exciting area of medical 
physics. However, EIT has not yet made the transition from an exciting medical 
physics discipline into wide spread routine clinical use. The technique still needs to 
break into widespread clinical acceptance and effort is continuing actively into the 
clinical trials and pilot studies which will achieve this. Technical advantages may 
allow us to obtain more accurate tissue characterisation and image quality and this 
will undoubtedly help to advance clinical acceptance. 

Outline of the thesis 

As already stated, the first part of the thesis is an introduction to the basic concepts 
of inverse and ill-posed problems. The aim of the thesis is not to study inverse 
problems in general so that this part is just a general introduction collected from 
literature [BaSTl IBB981 [Be86l \PS92i [Logg] . In the first chapter one can find the 
definitions of inverse and ill-posed problems, together with a few examples. Chapter 
2 presents some regularisation methods of ill-posed problems. First, the generalised 
solution of inverse problems is presented and then Tikhonov's regularisation method 
and the truncated singular value decomposition are considered. 

The second part of the thesis is concerned with an inverse problem in QCD: 
the determination of QCD condensates from r-decay data. Here, we will start in 
Chapter 3 with a theoretical description of r-decays: leptonic and hadronic decay 
width, hadronic polarisation function and its operator product expansion, and the 
dispersion relation satisfied by the polarisation function. This chapter represents 
the theory underlying the method we'll use to extract values for the condensates. 
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The experimental description of r-decays is disscusscd in Chapter 4. In Chapter 5 
we present a functional approach which allows us to extract within rather general 
assumptions values for the condensates from a comparison of the time-like r-dccay 
experimental data measured by the ALEPH Collaboration at LEP, with the asymp- 
totic space-like QCD prediction. Results for condensates of dimesion d = 4, 6, 8 and 
correlations between them for the V — A and V, A,V + A channels are presented in 
Chapters 6 and 7 respectively. 

The third part of the thesis is concerned with the inverse conductivity problem. 
In Chapter 8 we present the technique of EIT and its practical applications in 
medicine, industry and geophysics as well as a history of the problem. In Chapter 9 
we formulate the forward problem and aim to solve it by means of finite elements. 
Two reconstruction algorithms are discussed in Chapter 10. The first is based on 
reformulating the problem in terms of integral equations and aim to perform the 
reconstruction from a single set of measurements, while the latter is a linearisation 
type of algorithm which uses more sets of measurements. The second algorithm was 
also used on real data in Chapter 11, where reconstructions of a phantom immersed 
in a test tank filled with a conducting liquid are presented. Also measurements on 
the chest of a human object were taken and reconstructions performed with the aim 
of monitoring the lungs. 

There are also 4 appendices which contain supplementary mathematical infor- 
mations. In Appendix A we present the basic concepts of the theory of linear 
integral equations. In Appendix B we give the definition of Green's functions and 
we illustrate how they can be used to reduce the differential equations to integral 
ones. Appendix C presents the singular value decomposition method for integral 
operators. And last, in Appendix D, some basic concepts of statistics are presented: 
least-squares parameter estimation, confidence regions and the test of goodness- 
of-fit. 



Inverse problems 



Chapter 1 

Inverse and ill-posed problems 



Inverse problems of mathematical physics may be broadly described as problems of 
determining the internal structure or past state of a system from indirect measure- 
ments. Such problems would include for example the determination of diffusivities, 
conductivities, densities, sources, geometry of scatterers and absorbers and prior 
temperature distributions, to name just a few typical applications. Only recently a 
systematic treatment of such problems begun to emerge. The past few decades have 
witnessed a remarkable growth in inverse problems |5a00] . 

Inverse problems most often do not fulfil Hadamard's principle of well-posedness: 
they might not have a solution in the strict sense, solutions might not be unique 
and/or might not depend continuously on data. Hence their mathematical analysis 
is subtle. The belief of Hadamard that problems motivated by physical reality 
should be well-posed is essentially generated by physics of the nineteenth century. 
The requirements of existence, uniqueness and continuity of the solution are deeply 
inherent in the idea of a unique, complete and stable determination of the physical 
events. As a consequence of this point of view, ill-posed problems were considered, 
for many years, as mathematical anomalies and were not seriously investigated. 
The discovery of the ill-posedness of inverse problems has completely modified this 
conception. 

1.1 Inverse problems 

Suppose that we have a mathematical model of a physical process. We assume that 
this model gives a description of the system behind the process and its operating 
conditions and explains the principal quantities of the model (see Fig jl.ip : input, 
system parameters, output. 

In most cases the description of the system is given in terms of a set of equations 
(ordinary and/or partial differential equations, integral equations,...), containing 
certain parameters. 

The analysis of a given physical process via the mathematical model may be 
separated into three distinct types of problems. 

(A) The direct problem: Given the input and the system parameters, find out the 
output of the model. 

(B) The reconstruction problem: Given the system parameters and the output, 
find out which input has led to this output. 
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1. Inverse and ill-posed problems 



process 



input 


system 


output 




(parameters) 





Figure 1.1: Mathematical model of a physical process. 

(C) The identification problem: Given the input and the output, determine the 
system parameters which are in agreement with the relation between input 
and output. 

We call a problem of type (A) a direct (or forward) problem since it is oriented 
along a cause-effect sequence. In this sense problems of type (B) and (C) are called 
inverse problems because they are problems of finding out unknown causes of known 
consequences. It is immediately clear that the solution of one of the problems above 
involves a treatment also of the other problems. 

We give a mathematical description of the input, the output and the system in 
functional analytic terms. 

X space of input quantities; 

y space of output quantities; 

V space of system parameters; 

A{p) operator from X into y associated to p eV. 

In these terms we may reformulate the problems above in the following way: 

(A) Given / G A" and peV, find g = A{p)f. 

(B) Given g E y and p E V, solve the equation 

Af = g ifeX) (1.1) 
where A — A(p). 

(C) Given gey and f e X, find p eV such that 

A(p)f = g. (1.2) 

At first glance, the direct problem seems to be solved much more easier than the 
inverse problems. However, for the computation oi g = A{p)f it may be necessary to 
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solve a differential or integral equation, a task which may be of the same complexity 
as the solution of the equations in the inverse problem. 

In certain simple examples inverse problems can be converted formally into a 
direct problem. For example, if A has a known inverse then the reconstruction 
problem is solved by / = A~^g. However, the explicit determination of the inverse 
does not help if the output y is not in the domain of definition of A~^. This situation 
is typical in applications due to the fact that the output may be only imprecisely 
known and/or distorted by noise. 

As we stated above, a direct problem is a problem oriented along a cause-effect 
sequence; it is also very often a problem directed towards a loss of information: its 
solution defines a transition from a physical quantity with a certain information 
content to another quantity with smaller information content. In general it implies 
that the solution is much smoother than the data: the image provided by a band- 
limited system is smoother than the corresponding object, the scattered wave due 
to an obstacle is smooth even if the obstacle is rough, and so on. A more rigorous 
description of the loss of information typical for direct problems can be found in 
Section II. 3[ 

The conceptual difficulty common to most inverse problems is that by solving 
these problems, we would like to accomplish a transformation which should corre- 
spond to a gain of information. This provides the explanation of a typical mathe- 
matical property of inverse problems which is known as ill-posedness. 

1.2 Some examples of inverse problems 

Inverse problems fall mainly into three different but intimately related categories: 

• Inverse scaterring problems, 

• Inverse boundary value problems, 

• Inverse spectral problems. 

In the early 90's, the active research carried out in the field of inverse problems has 
brought a lot of new insight into the deeper nature of these problems and especially 
to the interrelation between them. In the following, we describe and discuss briefly 
the main features of each of these classes, give examples and further references. 

Inverse scattering problems 

Inverse scattering problems form undoubtedly one of the most studied set of inverse 
problems. The setting is the following: Far away from the target having unknown 
physical properties, a wave field is sent in. It is assumed that the interaction mech- 
anism of the wave field with the target is qualitatively known. The scattered field 
is measured, and from this data one attempts to reconstruct the properties of the 
scatterer. 



22 



1. Inverse and ill-posed problems 



A classical example of inverse scattering problems arises in quantum mechanics. 
Assume that we have a scattering potential q in M^. The quantum mechanical 
scattering with fixed energy E' = A;^, A; > = c = 1), is described by the 
Schrodinger equation 

{-/\-k'^ + q{x))i){x) = Q. (1.3) 

The potential q should decrease fast enough as x tends to infinity. The typical 
assumption about the field ip is that it is a superposition of the incoming plane 
wave and the scattered radiation field satisfying Sommerfeld's radiation condition 
at infinity, i.e., 

7/^(0;) = e^'=«-'^ + 7/;,c(a^), (1.4) 

where 

lim \x\{x ■ 'V — ik;)il)sc{x) = (1.5) 

|a;|^oo 

An equivalent way of formulating the radiation condition (11.51) is to assume that 

ik\x\ / 1 \ 

^sc(^) = ^A{x, e,k) + 0[^]. (1.6) 
\x\ \\^\ J 

The function A[x, 6, k) is called the scattering amplitude, and it is related to the 
scattering potential and the scattered field through 

A{x, e, k) = -^ [ e'''''-yq{yMy)d'y- (1.7) 

Depending on the type of measurements, one can now pose different inverse scatte- 
ring problems, of which we list the following: 

• Reconstruct the potential from the knowledge of the scattering amplitude at 
any energy 

{A{x,0,k) \ x,6 e S^, k eR+}; (1.8) 

• Reconstruct the potential from the knowledge of the scattering amplitude at 
fixed energy 

{A{x, 0,k) \x,0 e S^}, k>0 fixed; (1.9) 

• Reconstruct the potential from the knowledge of the backscattering amplitude 
at any energy 

{A{-e, e,k)\ees\ ke m+}. (i.io) 
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The first one of the above inverse problems is the most classical one. It is formally 
over-determined in the sense that the data set is indexed over a five dimensional 
space S'^ X S"^ X M_|_, while the unknown function q is over a three-dimensional space. 
Based on this over-determinacy, there is a rather simple way of seeing the uniqueness 
of the solution to this problem. Indeed, one can show that the scattering solution 
ip to (11. Sp behaves as 

ij{x) = e'''^-'' + O ( y] , asfc^oo. (1.11) 



Therefore, if we choose ^ G M'^ and let k tend to infinity while keeping the vector 
k{6 — x) = ^ fixed, we find that 

hm A{x,0,k) = J^[qm, (1.12) 

fc — *oo 

i.e., the scattering amplitude tends towards the Fourier transform of the potential. 
Therefore, the data of the problem determine the potential q uniquely. 

For the second problem posted, the over-determinacy of the data is one dimension 
less and consequently the problem to show the uniqueness of the solution is more 
difficult. As for the third one (called the inverse backscattering problem) we may 
mention that there are still open problems related to the uniqueness of the solution 
and reconstruction of the potential. 

Another classical, very important and closely related type of inverse scattering 
problems deals with obstacle scattering. A typical inverse obstacle scattering can 
be formulated as follows: Assume that in M^, there is an obstacle Q, whose shape 
one tries to recover from far field measurements. Assume that the medium outside 
the obstacle is governed by the equations of linear acoustics, i.e., the pressure field 
u satisfies the Helmholtz equation 

Au{x) + k^u{x) = 0, xeM.yn. (1.13) 

The pressure field is assumed to satisfy a boundary condition at dfl, typically the 
Dirichlet ( "soft sound" ) , Neumann ( "hard sound" ) or a mixed ( "impedance" ) condi- 
tion. For penetrable obstacles, the appropriate boundary condition is a transmission 
condition. Again, one probes the target by sending in an initial field uq, and the 
interacting total field is 

u{x) = uo{x) +Usc{x), (1.14) 

where the scattered field satisfies the outgoing radiation condition (II. 5p . or 

^ik\x\ / 1 \ 

u,,ix) = ——u^ix) + O { ] , (1.15) 
\x\ \|a;|^/ 

the function Uoo being the far pattern of the field which corresponds to the scattering 
amplitude here. The inverse scattering problem now consists in reconstructing the 
shape of the object from the far field patterns generated by a set of incoming fields. 
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Besides the acoustic inverse obstacle scattering, one can look at the similar prob- 
lem when the unknown target is illuminated with electromagnetic radiation. Since 
the electromagnetic fields satisfy the Helmholtz equation in vacuum, the transition 
from acoustics to electromagnetism does not seem so large. The boundary condi- 
tions, however, have vectorial nature and the field components will be coupled. 



Inverse boundary value problems 

Let us move now to the second large area of inverse problems, the inverse boundary 
value problems. The change in the setting when moving from inverse scattering 
problems to inverse boundary value problems is by no means sharp. Again, one has 
an object with unknown physical parameters and the objective is to find out these 
parameters in a non-invasive way. 

Let us start with a concrete example of impedance tomography: We ask whether 
it is possible to make an image of the internal electromagnetic structure of a body 
(e.g. human body) by injecting electric currents into the body and measuring the 
voltages needed to maintain the current. In practice, one attaches a number of 
electrodes on the surface of the body and measures the voltages needed to maintain 
the current configuration. By the linearity of the governing equations, dependence 
of the voltages on the currents is linear, so effectively the boundary data consists of 
a linear boundary map (or a matrix for the discretized version of the problem). 

The governing equation for the potential u in the body fl is simply the equation 
of continuity for the current j = a'Vu: 

V • aVu{x) = 0, xen, (L16) 

where (t{x) > is the conductivity distribution, assumed to be a scalar function of 
X. The current density through the boundary of the body is 



j{x) = a{x)^{x) 



(L17) 

an 



Assuming that the current density j is specified, one can solve the Neumann problem 
(11.161 - [T.17I) . In this way, one gets a complete collection of pairs 

{(j, M lac) I n satisfies (Ini-[LI7D}, (1.18) 

or, equivalently, one knows the Neumann-to-Dirichlet boundary map 



A 

dn 



^u\an. (1.19) 

an 



The inverse problem is to reconstruct a ixiVt from the knowledge of A. 

The close connection with the inverse scattering problems becomes more obvious 
if we modify Eq. (ll.l6p slightly. Let us introduce the function if) defined as 

i){x)=a{xf'^u{x). (1.20) 
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One can show that Eg. (11. 161) can be rewritten for ip as 

(A - q)%l) = 0, where q = ^ . (1.21) 

Thus, we are back to the Schrodinger equation with zero energy. 

As the discussion shows, there is a strong interrelation between the inverse sca- 
ttering and inverse boundary value problems. In many cases they can be shown to 
be equivalent: The knowledge of the boundary map determines the far field data 
uniquely and vice versa. 

Inverse boundary value problems get considerably more complicated if one allows 
for anisotropics in the medium. In fact, there are known limitations on the unique- 
ness for anisotropic inverse problems while the corresponding isotropic problems 
allow a unique solution. As an example, consider the anisotropic counterpart of 
Eq.dnS]), 

y —a'''—u = 0, inQe R", (1.22) 

^-^ OXi OXj 

l,J=l ■' 



the potential u satisfying the boundary condition 

dxj 

These equations can be written in coordinate free form using differential forms as 



= g\ (1.23) 

an 



dadu = inVL, (1.24) 



{adu)\Q^=g. (1.25) 

From this formulation, it can be shown that one can transform a by a diffeomorphism 
that leaves the boundary dVL untouched without affecting the boundary data. The 
question whether this is the only limitation on uniqueness is still open. 

Inverse spectral problems 

In the inverse spectral problems the input is the spectrum of an operator and one 
wishes to determine an unknown parameter of the operator. A classical example, 
also one of the simplest inverse problems in pure mathematics, is the one formulated 
by Mark Kac: "Can one hear the shape of the drum?" |Ka66] . Mathematically, the 
question is formulated as follows: Let be a simply connected, plane domain (the 
drumhead) bounded by a smooth curve 7, and consider the wave equation for u{x,t) 
(the displacement of the drumhead) on Q with a Dirichlet boundary condition on 7 
(the drumhead is clamped at the boundary): 

Au{x,t) = ———{x,t) mil, 

c (1.26) 

u{x, t) = on 7. 
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Looking for solutions of the form u{x,t) = Ree*'^*f(x) (normal modes) leads to an 
eigenvalue problem for the Dirichlet Laplacian on Q: 

Av{x) + Xv{x) = inQ, . 
v{x) = on 7, 

where A = u'^/c^. Kac's question means the following: is it possible to distinguish 
"drums" Qi and Q2 with distinct bounding curves 71 and 72, simply by "hearing" 
all of the eigenvalues of the Dirichlet Laplacian? 

Kac has showed that the asymptotic behaviour of (the resonance frequencies) 
at large k yields the volume and the total scalar curvature of Q or the length of 7 
|Ka66] . This kind of inverse problem has not been given real data applications, but, 
in some way, it is the pillars of the so-called geometric scattering theory. 



1.3 Ill-posed problems 

In the previous section we mentioned that a typical property of inverse problems is 
ill-posedness, a property which is opposite to that of well-posedness. 

The basic concept of a well-posed problem was introduced by the French math- 
ematician Jacques Hadamard in a paper published in 1902 on boundary-value pro- 
blems for partial differential equations and their physical interpretation |Hal902] . In 
this first formulation, a problem is called well-posed when its solution is unique and 
exists for arbitrary data. In subsequent work Hadamard emphasises the requirement 
of continuous dependence of the solution on the data |Ha23j , claiming that a solution 
which varies considerably for a small variation of the data is not really a solution 
in the physical sense. Indeed, since the physical data are never known exactly, this 
should imply that the solution is not known at all. 

From an analysis of several cases Hadamard concludes that only problems mo- 
tivated by physical reality are well-posed. An example is provided by the initial 
value problem for the D'Alembert equation which is fundamental in the description 
of wave propagation 

___(,.,)=0. (1.28) 

where c is the wave velocity. If we consider, for instance, the following Cauchy initial 
data at t = 

du 

m(x,0) = /(x), — (x,0) = 0, (L29) 

then there exists a unique solution given by 

u{x, t) = ^ [fix - at) + fix + at)] . (L30) 

This is a solution for any continuous function fix). Moreover it is obvious that a 
small variation of fix) produces a small variation of uix,t). 
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The previous problem is well-posed and, of course, basic in the description of 
physical phenomena. It is an example of a direct problem. An impressive example 
of an ill-posed problem and, in particular, of a non- continuous dependence on the 
data, was also provided by Hadamard |Ha23] . This problem which, at that time, 
was deprived of a physical motivation, is the Laplace equation in two variables 

+ 1^(^,2/) =0. (1.31) 

If we consider the following Cauchy initial conditions at ?/ = 

M a;,0 = -COS nx , — a;, = 0, 1.32 
n ay 

then the unique solution is given by 

u{x,y) = — cos{nx) cosh{ny). (1.33) 
n 

The factor cos(nx) produces an oscillation of the surface representing the solution of 
the problem. This oscillation is imperceptible near y = but becomes enormous at 
any given finite distance from the x-axis when n is sufficiently large. More precisely, 
when n — oo, the data of the problem tend to zero but, for any finite value of y, 
the solution tends to infinity. 

This is now a classical example illustrating the effects produced by a non- 
continuous dependence of the solution on the data. If the oscillating function de- 
scribes the experimental errors affecting the data of the problem then the error 
propagation from the data to the solution is described by Eq. fll.33p and its effect 
is so dramatic that the solution corresponding to real data is deprived of physical 
meaning. Moreover it is also possible to show that the solution does not exist for 
arbitrary data but only for data with specific analyticity properties. 

A problem satisfying the requirements of existence, uniqueness and continuity is 
now called well-posed in the sense of Hadamard, even if the complete formulation 
in terms of the three requirements was first given by R. Courant |CH53b] . The 
problems which are not well-posed are called ill-posed or also incorrectly posed or 
improperly posed. Therefore an ill-posed problem is a problem whose solution is not 
unique or does not exist for arbitrary data or does not depend continuously on the 
data. 

The previous observations and considerations can justify now the following ge- 
neral statement: A direct problem, i.e., a problem oriented along a cause-effect 
sequence, is well-posed while the corresponding inverse problem, which implies a re- 
versal of the cause-effect sequence, is in general ill-posed. This statement, however, 
is meaningful only if we provide a suitable mathematical setting for the description 
of direct and inverse problems. 

The first point is to define the class of objects to be imaged, which will be 
described by suitable functions with certain properties. In this class we also need 
a distance, in order to establish when two objects are close and when they are not. 



28 



1. Inverse and ill-posed problems 



In such a way our class of objects takes the structure of a metric space of functions. 
We denote this space by X and we call it the object space. 

The second point is to solve the direct problem, i.e. to compute, for each object, 
the corresponding image which can be called the computed image or the noise-free 
image. Since the direct problem is well-posed, to each object we associate one, and 
only one, image. As we already mentioned, this image may be rather smooth as a 
consequence of the fact that its information content is smaller than the information 
content of the corresponding object. This property of smoothness, however, may not 
be true for the measured images, also called noisy images, because they correspond 
to some noise-free image corrupted by the noise affecting the measurement process. 

Therefore the third point is to define the class of the images in such a way that 
it contains both the noise-free and the noisy images. It is convenient to introduce 
a distance also in this class. We denote the corresponding function space by y and 
we call it the image space. 

objects with very 




same image 

Figure 1.2: Schematic representation of the relationship between objects and 
images. 

In conclusion, the solution of the direct problem defines a mapping (operator), 
denoted by A, which transforms any object of the space X into a noise-free image 
of the space 3^. This operator is continuous, i.e. the images of two close objects are 
also close, because the direct problem is well-posed. The set of the noise-free images 
is usually called, in mathematics, the range of the operator A, and, as follows from 
our previous remark, this range does not coincide with the image space y because 
this space contains also the noisy images. 

By means of this mathematical scheme it is possible to describe the loss of 
information which, as we said, is typical in the solution of the direct problem. It 
has two consequences. First, it may be possible that two, or even more, objects 
have exactly the same image. In the case of a linear operator this is related to 
the existence of objects whose image is exactly zero. These objects will be called 
invisible objects. Then, given any object of the space X, if we add to it an invisible 
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object, we obtain a new object which has exactly the same image. Secondly, and 
this fact is much more general than the previous one, it may be possible that two 
very distant objects have images which are very close. In other words there exist 
very broad sets of distinct objects such that the corresponding sets of images are 
very small. All these properties are illustrated in Fig J1.2[ 

If we consider now the inverse problem, i.e. the problem of determining the 
objects corresponding to a given image, we find that this problem is ill-posed as 
a consequence of the loss of information intrinsic to the solution of the direct one. 
Indeed, if we have an image corresponding to two distinct objects, the solution of 
the inverse problem is not unique. If we have a noisy image, which is not in the 
range of the operator A, then the solution of the inverse problem does not exist. 
If we have two neighbouring images such that the corresponding objects are very 
distant, then the solution of the inverse problem does not depend continuously on 
the data. 

1.4 A few examples of ill-posed problems 

Predholm integral equations of the first kind 

A Fredholm integral equation of the first kind is an equation of the form 



where y is a given function (usually called the data), K{-,-) is the kernel of the 
equation and the solution x is the unknown function which is sought. There are 
several observations concerning this equation. The first one is that the function y 
inherits some of the smoothness of the kernel K and therefore a solution may not 
exist if y is too roughly behaved. For example, if the kernel K is continuous and x 
is integrable, then the function y defined by Eq. fl 1.341) is also continuous and hence 
if the given function y is not continuous while the kernel is, then Eq. fl 1.34p can not 
have an integrable solution. Consequently, the question of existence of solutions is 
not trivial and requires more detailed knowledge of the properties of K. 

Another point to be considered is the uniqueness of solutions. For example, if 
K{s,t) = ssint, then the function x{t) = 1/2 is a solution of 



but so is each of the functions x„(t) = 1/2 + sin(nt), for n = 1, 2, 3, .... 

A more serious concern arises from the Riemann-Lebesgue lemma which states 
that if K{-, ■) is any square integrable kernel, then 
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From this it follows that if x is a solution of Eq. fll.34p and A is arbitrary, then 

K{s,t) {x(t) + Asm{nt)) dt y{s) as n — *• oo. (1-37) 

Therefore for large values of n the slightly perturbed data 

y{s)=y{s) + A[ K{s,t)sm{nt)dt (1.38) 
Jo 

corresponds to a solution x{t) + Asin{nt) which differs markedly from x{t). Hence, 
for Fredholm equations of the first kind, solutions generally depend discontinuously 
upon the data|l| 

Cauchy problem for the Laplace equation 

The simplest example of ill-posed problems for the Laplace equation is a mixed 
boundary value problem in two dimensions. The problem is to determine in a 
rectangle {0<x<7r, < y < yo} a function of two variables u{x, y) satisfying the 
following conditions 



Au{x,y) = 0, 

u{0,y) = u{7T,y) = 0, , . 

«(x,0) = /o(x), ^'-^^^ 

Uy{X,0) = fl{x). 

The solution of the Laplace equation satisfying the homogeneous conditions on the 
edges of the strip can be represented in the form 

oo 

u{x, y) = J2 K^*'^ + ^kC'^^) sin kx. (1.40) 
fc=i 

From the initial conditions one finds 



TV 



~ (^J foix) sinkxdx + — J f i(x) sin kxdx^ , (1-41) 



bk = — \ fo{x)sinkxdx — 7 / fi{x)smkxdx ) . (1-42) 



^ \Jo k Jo 

Thus, the solution of the problem (11.391) is unique but its existence is not guaranteed 
and depends on the integrability properties of /o and /i and on the properties which 
u{x, y) should have. 

Let us now consider the following solutions of f ll.39p 

Un{x,y) = a„e"^ sin rax. (1-43) 



^For more details related to the solutions of the Fredholm integral equations of the first kind see 
Appendix A. 
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It is clear that the Cauchy data which lead to these solutions are: 

fon = anSmnx, (1-44) 

fin = nan sin nx. (1-45) 

Obviously, an appropriate choice of n and a„ may render these Cauchy data arbi- 
trarily small in the norn|§, while the function Un will be arbitrarily large for any 
fixed y. 

Analytic continuation 

There are several settings of the analytic continuation problem. We present here 
one classical example of analytic continuation for functions of a complex variable. 

Let f{z) be an analytic^] function of a complex variable, which is reg ulai0 within 
a bounded region Q on the complex plane and continuous in the closure Q, 

\f{z)\<c, zen. (1.46) 

Let dQ be the boundary of Q, Li, La be parts of dQ: Li U Fa = dQ, Fi n Fa = 0. 
Suppose also that f{z) is specified on Fi, and the problem is to determine f{z) 
within the interior part of Q. If dQ = Fi, the solution is provided by the Cauchy 
integral. 

If Fi 7^ dQ, then to find the analytic continuation of f{z) is equivalent to the Cauchy 
problem for the Laplace equation. 

Let the value of a harmonic function u and its normal derivative du/dn be given 
on Fi. Denote by f{z) the function 

f{z)=u{z)+iv{z), z = x + iy, (1-48) 

where v is the function conjugate to u. It is known that on Fi 

v{z) = [ ■^u{s)ds + Cu (1.49) 
Jz^ on 



where zq is one of the ends of Fi and Ci a constant. Hence, if u{z), du{z)/dn are 
known on Fi, then the analytic function f{z) may be deemed given on Fi. 



^Independent of the specific choice of the norm. 

■^An analytic function is an infinitely differentiable function such that the Taylor series at any point 



xq in its domain, 



n=0 



is convergent for x close enough to xq and its value equals /(x). 

''a function is termed regular if and only if it is analytic and single-valued throughout a region fi. 
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From the Caucliy-Riemann conditions it follows that 



du 
dn 



dv 
ds 



;i.5o) 



where dv/ds is the derivative of v along Fi. Taking the derivative of f{z) and f{z) 
along Fi yield^ 



du 
dn 



ld_ 
2ds 



if - f) 



X51] 



and we arrive at the Cauchy initial conditions for u{x,y). 

Thus, if Fi 7^ dfl, the problem of analytic continuation is equivalent to the 
Cauchy problem for the Laplace equation and so is ill-posed. 



1.5 How to cure ill-posedness 

The property of a non-continuous dependence of the solution on the data strictly 
applies only to ill-posed problems formulated in infinite dimensional spaces like the 
ones discussed in the previous section. In practice one has discrete data and one has 
to solve discrete problems. These, however, are obtained by discretizing problems 
with very bad mathematical properties. What happens in these cases? 

If we consider a linear inverse problem, its discrete version is a linear algebraic 
system, apparently a rather simple mathematical problem. Many methods exist for 
solving numerically this problem. However the solution often does not work. A 
description of the first attempts of data inversions is given by S. Twomey in the 
preface of his book |Tw77] : 'The crux of the difficulty was that numerical inversions 
were producing results which were physically unacceptable but were mathematically 
acceptable (in the sense that had they existed they should have given measured 
values identical or almost identical with what was measured)'. These results were 
'rejected as impossible or ridiculous by the recipient of the computer's answer. And 
yet the computer was often blamed, even though it had done all that had been 
asked of it'. ...'Were it possible for computers to have ulcers or neuroses there is 
little doubt that most of those with which early numerical inversion attempts were 
made would have required both affiictions' |Tw77j . 

The explanation can be found having in mind the examples discussed in the 
previous section, where small oscillating data produce large oscillating solutions. 
In any inverse problem, data are always affected by noise which can be viewed 
as a small randomly oscillating function. Therefore the solution method amplifies 
the noise producing a large and wildly oscillating function which completely hides 
the physical solution corresponding to the noise-free data. This property holds 
true also for the discrete version of the ill-posed problem. Then one says that the 
corresponding linear algebraic system is ill-conditioned: even if the solution exists 



^Here, / denotes the complex conjugate of /. 
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and is unique, it may be, and is in general, completely corrupted by a small error 
on the data. 

In conclusion, we have the following situation: we can compute one, and only 
one, solution of our algebraic system but this solution may be unacceptable for the 
reasons indicated above; the physically acceptable solution we are looking for is not 
a solution of the problem but only an approximate solution in the sense that it does 
reproduce the data not exactly but only within the experimental errors. However, 
if we look for approximate solutions, we find that they constitute a set which is 
extremely broad and contains completely different functions, a consequence of the 
loss of information in the direct problem. Then the question arises: how can we 
choose the good ones? 

We can state now the 'golden rule' for solving inverse problems which are ill- 
posed: search for approximate solutions satisfying additional constraints coming 
from the physics of the problem. 

The set of the approximate solutions corresponding to the same data function is 
just the set of objects with images close to the measured one. The set of objects is 
too broad, as a consequence of the loss of information due to the imaging process. 
Therefore we need some additional information to compensate this loss. This infor- 
mation, which is also called a priori or prior information, is additional in the sense 
that it cannot be derived from the image or from the proprieties of the mapping A 
which describes the imaging process but expresses some expected physical proper- 
ties of the object. Its role is to reduce the set of the objects compatible with the 
given image or also to discriminate between interesting objects and spurious objects, 
generated by uncontrolled propagation of the noise affecting the image. 

The idea of using prescribed bounds to produce approximate and stable solutions 
was introduced by C. Pucci in the case of the Cauchy problem for the Laplace equa- 
tion jPu55] . i.e., the first example of an ill-posed problem discussed by Hadamard. A 
general version of similar ideas was formulated independently by V.K. Ivanov |lv62] . 
His method and the method of D.L. Phillips for Fredholm integral equations of the 
first kind |Ph62] were the first examples of regularisation methods for the solution of 
ill-posed problems. The theory of these methods was formulated by A.N. Tikhonov 
one year later [Ti63] . 

The principle of the regularisation methods is to use the additional information 
explicitly, at the start, to construct families of approximate solutions, i.e. of objects 
compatible with the given image. These methods are now one of the most powerful 
tools for the solution of inverse problems, another one being provided by the so-called 
Bayesian methods, where the additional information used is of statistical nature. 

We will continue the discussion on regularisation methods in the next chapter. 
First, we will present the generalised solution of inverse problems, describe the reg- 
ularisation method introduced by Tikhonov and also the truncated singular value 
decomposition. In the last section, we will give a general definition of a regularisa- 
tion algorithm and give some properties a regulariser should have so that it gives 
approximate and stable solutions to inverse problems. 



Chapter 2 

Regularisation of ill-posed 

problems 



2.1 The generalised solution 

Given the noisy image g & y and the hnear operator A describing the imaging 
system, we are interested in solving the hnear equation 

Af = g (2.1) 
for f ex. 

We assume that the operator A has a singular value decomposition so that we 
can write (see Appendix [C]) 

oo 

Af = Y,^Af,^j)^^j, (2.2) 

i=i 

where aj are the singular values of A and vj, uj are singular functions on X and y 
respectively. 

The problem (12. ip is, in general, ill-posed in the sense that the solution is not 
unique, does not exist, or else, does not depend continuously on the data. 

Uniqueness does not hold when the null-space of the operator A, J^{A), i.e. the 
set of the invisible objects / such that Af = 0, is not trivial. 

The procedure most frequently used for restoring uniqueness is the following one. 
Any element / of the object space X can be represented by 

oo 

f = Y.^f.Vj)xv, + v, (2.3) 
i=i 

where v is the projection of / onto A/'(A), while the first term is the component of / 
orthogonal to M{A). The term v can be called the invisible component of the object 
/ because it does not contribute to the image of /, Af . Since the invisible component 
cannot be determined from Eq. fl2.1l) . it may be natural to look for a solution of this 
equation whose invisible component is zero. Such a solution is unique because from 
Eqs. (12.21) and (12. 3p . with t> = 0, we easily deduce that Af = implies / = 0. If this 
solution exists, it is denoted by /"*" and called minimal norm solution. Indeed, any 
solution of fl2.ip is given by 

f = f'- + v, (2.4) 
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where v is an arbitrary element of //{A). Since v is orthogonal to we have 

ll/ll^ = lirilUll^ll^ (2.5) 

and therefore the solution with f = 0, i.e., /+, is the solution of minimal norm. 

As concerns the existence of a solution of (12. ip and, in particular, of we first 
have to distinguish between the following two cases. 

• The null space of the adjoint A* of A, //{A*), contains only the zero element. 
Then the singular functions (vectors) Uj constitute an orthonormal basis in 3^ 
and any noisy image g can be represented by 

oo 

g = ^{g,Uj)yUj. (2.6) 
i=i 

By comparing this representation with the SVD of A, Eg. (12.21) . we see that a 
solution of (12. ip may exist. 

• The null space ^/{A*) contains non-zero elements. In such a case the singular 
functions (vectors) Uj do not constitute an orthonormal basis in y and the 
noisy image g can be represented as follows 

oo 

g = ^{g,Uj)yUj + u (2.7) 
i=i 

where u is the component of g in Af{A*), i.e. the component of g orthogonal to 
the range of A. Notice that, if the mathematical model of the imaging system 
is physically correct, the presence of this term is an effect due to the noise. If 
M 7^ 0, by comparing the representation (12.20 of Af with the representation 
(12. 7p of g, we see that there does not exist any object / such that Af coincides 
with g. Then we can look for objects / such that Af is as close as possible to 
/, i.e. for objects which minimise the discrepancy functional 

1 1^/ ^ S'l |y = '^i^i^um. (2.8) 

Any solution of this variational problem is called a least-squares solution. 

The concept of least-squares solutions is more general than the concept of solu- 
tion because a solution of (12.10 is also a least-squares solution. More precisely, the 
set of the least-squares solutions coincides with the set of the solutions if and only 
if the minimum of the discrepancy functional (12. 8p is zero. This remark shows that, 
without loss of generality, we can investigate the problem of existence in the case of 
the least-squares solutions. 

Solving problem (12.80 is equivalent to solving its Euler equation, which is given 

by 



A*Af = A*g. 



(2.9) 



2.1. The generalised solution 



37 



From the SVD of the operator A, Eq. fl2.2p . and of the operator A* 

oo 

A*9 = Y.''^^9,u,)yv„ (2.10) 
i=i 

we obtain 

oo 

A*Af = J2^]U\vjUv,. (2.11) 

i=i 

If we insert these representations into Eg. 02.91) and compare the coefficients of Vk, 
we find that the components of any solution / of fl2.9p are given by 

Crj{f,Vj)x = (Tj{g,Uj)y, (2.12) 

and therefore 

{f,Vj)x = —{g,Uj)y. (2.13) 

In such a way the existence of least-squares solutions has been reduced to the 
existence of elements of the object space X whose components with respect to the 
singular functions Vj are given by Eq. fl2.13p . 

According to (12.130 we can introduce the following formal solution 

oo ^ 

r = J2-i9,^j)yVj. (2.14) 

We say that this solution is formal because it is given by a series expansion so that 
the solution exists if and only if the series is convergent. 

If we consider the convergence in the sense of the norm on X, then this con- 
vergence is assured if and only if the sum of the squares of the coefficients of the 
eigen-functions Vj is convergent. We obtain the following condition 



1 2 

■^\i9,Uj)y\ < oo (2.15) 



E 

j=i "J 

which is also called the Picard criterion for the existence of solutions or least-squares 
solutions of the linear inverse problem we are considering |Gr93] . 

It is important to point out that, if the singular values aj accumulate to zero 
then condition (12.151) may not be satisfied by any arbitrary noisy image g. If the 
condition (I2.15P is not satisfied, then no solution or least-squares solution of the 
inverse problem exists. 

The functions g satisfying the Picard criterion are images in the range of A, 
TZ{A). For any one of these functions, the series (12.140 defining /"*" is convergent. 
Then we can conclude that for any image g satisfying the Picard criterion there 
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exists a unique generalised solution f^, whose singular function expansion is given 
by Eq.^^. 

The generalised solution defines a generalised inverse operator . This operator 
is not defined everywhere on 3^ but only on the set of functions g satisfying the Picard 
criterion. This set is the domain of the operator , ^{A^). Moreover the operator 

is not continuous or, in other words, the generahsed solution /+ does not depend 
continuously on the image g. In order to prove this statement, let us assume that g 
is an image satisfying the Picard criterion and let us consider a sequence of images 
given by 

9j=9 + V^Uj. (2.16) 

If we assume that the singular values tend to zero, which is in general the case, it is 
obvious that 

\\9j - 9\\y = ^ 0- (2.17) 

On the other hand, if we denote by /"*" the generalised solution associated with g, 
and by fj~ the generalised solution associated with gj, we have 



f; = r + ^v, (2.18) 

so that 

11/^-/^11 = (2-19) 



In such a way we have found a sequence of images, converging to g, such that the 
sequence of the corresponding generalised solutions does not converge to /"*". 



2.2 Tikhonov's regularisation method 

The generalised solution of an ill-posed or ill-conditioned problem is not physically 
meaningful because it may be completely corrupted by the noise propagation from 
data to solution. For this reason we must look for approximate solutions satisfying 
additional constraints suggested by the physics of the problem and regularisation is 
a way for obtaining such solutions. 

The starting point is to define a family of regularised solutions fx, depending 
on a regularisation parameter A > 0, as the family of the functions minimising the 
functionals 

^,{f;g) = \\Af-g\\l + X\\f\\l (2.20) 

where g is the given image. The meaning of the regularisation parameter A will 
become clear in Section [231 where a general theory of regularisation algorithms will 
be presented. Let us first discuss how can one find a unique solution to the inverse 
problem (12.11) which simultaneously minimise the functional (I2.20p . 



2.3. Truncated SVD 
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The Euler equation associated with the minimisation of this functional is given 

by 

{A*A + XI)f = A*g. (2.21) 
An object is a minimum point fx of the functional (12.201) if and only if it is a solution 

of dai]). 

In order to solve this equation, let us represent an arbitrary element / of A* in 
terms of the singular functions Vj of the operator A and of the elements v (orthogonal 
to all Vj) of the null space of A, as in Eq. (l2.3p . If we insert this representation in 
(I2.2ip and we take into account Eqs. (12.101) and (12.111) . we obtain 

oo oo 

^^(aj + A)(/, Vj)xVj + Xv = J2 ^A9^ ^j)y^r (2-22) 
i=i i=i 

It follows that there exists a unique solution fx of (12.211) . which can be obtained 
from (12. 3p with v = and with coefficients (/, Vj)x given by 

{a] + X){f,v,)^ = a,{g,u,)y. (2.23) 

In conclusion we find 

oo 

h = J2^Tx(9,n,)yv,. (2.24) 

The series at the r.h.s. of this equation does always converge (thanks to the factors 
aj, the coefficients tend to zero more rapidly than the components {g, Uj)y of g) and 
therefore the regularised solution fx exists for any noisy image g. 

2.3 Truncated SVD 

The representation (12.240 of the regularised solution can be recast in the following 
form 

fx = j2^(9,u,)yvj, (2.25) 

where 



■ 1 



W., = (2.26) 

This expression shows that the regularised solution fx can be obtained by a filtering 
of the singular value decomposition of the generalised solution: the components 
of corresponding to singular values much larger than A are taken without any 
significant modification, whereas the components corresponding to singular values 
much smaller than A are essentially removed. 
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One possibility is to replace the smooth filter given in Eg. (12.261) by a sharp one, 
i.e., to take in the singular function expansion of the generalised solution only the 
terms corresponding to singular values greater than a certain threshold value. Since 
the singular values are ordered to form a decreasing sequence, those greater than 
the threshold value are those corresponding to values of the index less than a certain 
maximum integer. 

Let us denote by J the number of singular values satisfying the condition 
(t|>A, j<J, (2.27) 
then the approximate solution provided by the truncated SVD is as follows 
J 1 

Fj = Y.-{g^u^)yVj. (2.28) 

This equation can be obtained from (12.251) by taking Wxj = 1 when (Jj>X and 
Wxj = when crj < A. 



2.4 Regularisation algorithms 

The methods investigated in the previous sections can be embedded in a more general 
approach called by Tikhonov the regularisation method |Ti63i [TA77j . It consists in 
the introduction of families of continuous approximations to the generalised inverse 
of the operator A. 

A one parameter family of operators {-Ra}a>o is called a regularisation algorithm 
or a regulariser for the solution or general solution of Eq. Ol.ip . if the following 
conditions are satisfied |TA77] : 

i) for any A > 0, -Ra : 3^ — ^ is continuous; 

ii) for any g such that g G 7^(v4) 

lim||i?A^7-/+||;t = (2.29) 

where is the generalised solution (or solution when exists) of Eq. fll.ll) . 

When Rx is linear we have a linear regularisation algorithm; the variable A is 
called the regularisation parameter. 

Eq. (l2.29p can also be written as follows 

hm\\RxAf+ - f+\\;r = 0. (2.30) 
Therefore the operator Tx : X ^ X , defined by 



Ta = RxA 



(2.31) 



2.4. Regular isat ion algorithms 
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is an approximation of the orthogonal projection onto N'{A)-'- or of the identity 
operator when A^(^) = {0}. 

Conditions i) and ii) imply that, for any A and for any exact data g G TZ{A), 
R\g is a continuous approximation of the solution or generalised solution of Eg. (II .11) . 
However, the important case is that of noisy data Qe-, with ^ TZ{A), which are 
close to exact data g G Ti^A), — (^Hy < if e is small enough. In this case, no 
solution of the equation Af = may exist. 

If we consider the functions fx = Rxge: it is easy to see that there must exist 
an optimum value of A such that fx is as close as possible to f~^, the generalised 
solution associated with the exact data g. Assuming that Rx is linear and that the 
noisy data gs are written in the form g^^ = Af + Ws, we get 

Rxge - /+ = iRxAf+ - /+) + RxWe, (2.32) 

and therefore 

I \Rxge -f^\\x< ^(A; /+) + eN{X), (2.33) 

where 

iV(A) = ||i?A||, ^(A;/+) = ||i?A^/+-/+|U and e = \\w,\\y. (2.34) 

The first term of the r.h.s. in Eq. fl2.33p represents the approximation error intro- 
duced by the choice of a non-zero value of the regularisation parameter; it tends to 
zero when A — >■ 0. The second term represents the error on the approximate solu- 
tion induced by the error on the data; it tends to infinity when A ^ 0. Therefore it 
is necessary to find a compromise between approximation and error magnification. 
Assume, for simplicity that iV(A) and a;(A; /+) are monotonous functions of A (this 
condition is satisfied by all the regularisation algorithms used in practice) and, more 
precisely, that A^(A) is a decreasing function, with A^(O-I-) = +oo, while u;(A; /"*") is 
an increasing function, with uj{0+; f^) = 0. Under these assumptions, there exists 
a unique value of A, A(£:), which minimises the r.h.s. of Eq. (12.331) and which re- 
presents the optimum compromise between approximation and error magnification. 
Furthermore X = X{e) ^ 0, when £ — >• 0, and Rxge ~^ f^- 

The previous argument implies that a regularisation algorithm can give appro- 
ximate and stable solutions which converge to the exact solution when the error on 
the data tends to zero. 

The Tikhonov regulariser 

The most intensively investigated example of a regularisation algorithm is the so- 
called Tikhonov regulariser, given by 

Rx = {A*A + XI)-^A\ (2.35) 

Its remarkable properties derive from the fact that it can be obtained by minimising 
the functional fl2.20p (see Section [272]) . In particular it is easy to show that Rxg G 
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Af{A)-^, for any g E y and this implies that {-Ra}a>o is a regularisation algorithm 
for /+. Using the spectral representation of A* A, indeed, it is easy to show that 

u{X- /+) = \\RxAt -nU = All (AM + A/)-V+IU (2.36) 

tends to zero when A — 0. Furthermore, for any /"*", uj{\] /+) is an increasing 
function of A. 

It is also easy to find an estimate for iV(A), Eq. fl2.34p . 

\\Rxg\\l = {AR^g , {A*A + \I)-^g)y<^\\g\\l (2.37) 
since 1 1 AR\ 1 1 < 1 , and therefore 

N{\) < -L. (2.38) 

From this bound and Eq. fl2.33p one can derive results concerning possible choices of 
the function A(£:). 

Spectral windows 

The regulariser fl2.35p can be written in the following formal way 

= Ux{A*A)A*, (2.39) 

with the help of the function 

UM = {fi + X)-\ (2.40) 

This remark suggests a way for defining a wide class of regularisation algorithms 
which can be applied whenever the spectral representation of A* A is known |Ba65l 

Consider a family of real-valued, piecewise continuous functions {U\}xyo, defined 
on the interval [0, ||A|p] and assume that they satisfy the following conditions: 

i) for each A > 0, there exists a constant cx such that 

|f/A(/i)| <ca; V/iG [0,P||2]; (2.41) 

ii) for each A > 

< /xf/A(/i) < 1; y fie[0,\\A\\']; (2.42) 



iii) 



lim/if/A(/i) = l; V/xG (0,P||2]. (2.43) 

A \0 
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Then the family of operators Rx, defined by Eq. fl2.39p is a regularisation algorithm 
for /"•". This result derives from the following remarks. Thanks to condition i), for 
each A > 0, the operator Rx is bounded. Furthermore the following relation holds: 

UxiA*A)A* = A*UxiAA*) (2.44) 

since it is true for polynomials and therefore it is also true for continuous functions. 
This relation implies that, for any A > and any g & y 

fx = Rx9 e M{A)^. (2.45) 

Finally, from the spectral representation of A* A, from conditions ii) and iii) and the 
dominated convergence theorem^, property fl2.30p can be derived. 

The function Ux{fi) defined in Eq. fl2.40p satisfies the previous conditions. Ano- 
ther important example is given by 

"^(■^)-{% : r>f i^*) 

For a regularisation algorithm of the type (12.390 . the operator Tx, Eq. (l2.3ip . is 

Tx = Wx{A*A), (2.47) 

where WxifJ^) = fiUxifi)- The function Wx{fi) can be called a spectral window. The 
justification of this name derives from the example 02.461) since in such a case the 
function Wx{fi) is zero in a neighbourhood of and is unity elsewhere. 

In the case of a compact operator, the regulariser (12.391) can be expressed in 
terms of the singular value decomposition 

CXD 

Rxg = Y,f'k'^'^M{9,Uk)vk. (2.48) 

In particular, in the case (|2.46p we have 

Rxg = J2^k'^^(9,Uk)vk. (2.49) 

and this is the well-known method of the truncated singular value decomposition 
discussed in Section 12.31 



^The dominated convergence theorem states: Let X he a measurable space and let g, /i, /2, ... be 
measurable functions such that g < oo and |/„| < g for each n. If /« ^ / almost everywhere, then 
/ is integrable and 

lim / /„ = / /. 
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2.5 Choice of regularisation parameter 

The choice of the regularisation parameter is a crucial and difficult problem in the 
theory of regularisation. This point has been widely discussed in the mathemati- 
cal literature. No precise recipe has been discovered which could be used for any 
problem. The existence of a recipe depends a lot on the application and the actual 
information one has at hand. In what follows we will restrict ourselves to the case 
of Tikhonov's regularisation algorithm discussed in Section 12.21 and summarise the 
main methods which are used in practice. 

As we know from the discussion of Section 12. 4[ for any image g there exists 
an optimum value, Aopt, of the regularisation parameter. For this value of A, the 
corresponding regularised solution fx has minimal distance from the true object /. 
The problem is that the determination of this optimal value requires the knowledge 
of/. 



Regularised solution with prescribed energy 

If we do not know / but know its norm E = ||/||, also called 'energy', one may try 
the value A = \{E) such that the corresponding regularised solution has the same 
energy as the true object, i.e. 

\\fx\\=E. (2.50) 

X{E) is a decreasing function of E. Therefore, if we overestimate the energy of 
the object, we obtain a value of the regularisation parameter smaller than that 
corresponding to the exact energy of the object. In such a case the regularised 
solution will show a higher noise contamination. 



Regularised solution with prescribed discrepancy 

If we know a precise estimate e of the energy of the noise, then the estimate is the 
value A = A(£:) such that the discrepancy of the corresponding regularised solution 
is just equal to e, i.e. 

ll^/A-^?||=e. (2.51) 

This method is also known as the discrepancy principle |Mo84] . X{e) is an increasing 
function of e. Therefore, if we overestimate the energy of the noise, we get a value 
of the regularisation parameter which is larger than that corresponding to the exact 
energy of the noise. In such a case the regularised solution will show a smaller noise 
contamination. 



The Miller method 

An approach to ill-posed problems proposed by Miller |Mi70] can also be considered 
as a method for estimating the value of the regularisation parameter. In this appro- 
ach it is assumed that one knows both a bound on the energy and a bound on the 
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discrepancy of the unknown object /. Then the set of all objects / satisfying the 
two conditions 

\\Af-g\\^<e\ \\f\\^<E\ (2.52) 

is called the set of admissible approximate solutions. This set corresponds to the 
intersection of the ball of the objects with squared energy smaller than E"^ and of 
the ellipsoid of the objects with discrepancy smaller than e^. If this intersection is 
not empty, then the pair {e, E} is said to be permissible. 

Now, if the pair {e, E} is permissible, it has been shown by Miller |Mi70] that 
the regularised solution corresponding to the following value of the regularisation 
parameter 

A = (2.53) 

satisfies the condition fl2.52l) except for a factor of \/2 and therefore it is essentially 
an admissible approximate solution. 



Generalised cross-validation 

The methods considered previously require the knowledge of e or of or of both. 
In many cases one does not have a sufficiently accurate estimate of these quantities 
and therefore it is important to have methods which do not require this kind of 
information. One such method is that of cross-validation |Wa77j which can only 
be used in problems with discrete data. It is based on the idea of letting the data 
themselves choose the value of the regularisation parameter. In other words one 
requires that a good value of the regularisation parameter should predict missing 
data values. 

The mathematical formulation of this method is more complicated and will not 
be presented here. For more details one can look in |BB98] . 



L-curve method 

This graphically motivated method, introduced by Hansen |Ha92j . is another method 
which does not require information about the energy of the noise or of the true object. 
The starting point is the plot of E{fx) versus e{fx] g), introduced in connection with 
Miller's method. This curve has, in many cases, a rather characteristic L-shaped 
behaviour in a log- log plot. 

A qualitative explanation of this behaviour is the following. We recall that E{f\) 
is large for small A and small for large A while e{fx\ g) has the opposite behaviour. 
Therefore E{fx) is large when e{fx; g) is small, and conversely. This is the trade-off 
between noise-propagation error and approximation error. Now the vertical part of 
the L-curve corresponds to values of the regularisation parameter such that fx is 
dominated by the noise propagation error. As a consequence E{fx) is very sensitive 
to variations of A while e{fx', g) is not. Analogously the horizontal part of the L-curve 
corresponds to values of the regularisation parameter such that fx is dominated by 
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the approximation error. As a consequence e{fx; g) is very sensitive to variations of 
A while E{fx) is not. 

Now the L-curve method consists of taking as an estimate of the regularisation 
parameter the value of A corresponding to the corner of the L-curve. In fact this 
point should correspond to the compromise between approximation error and noise- 
propagation error. From the computational point of view a convenient definition of 
the L-curve corner is the point with maximum curvature. 

The L-curve method, even if it can be useful in some cases, does not work 
in all cases and has some theoretical and practical inconveniences. It has been 
shown |EG94| IVo96] that, in certain cases, it does not provide a regularised solution 
converging to the exact one when the noise tends to zero. Moreover, examples can 
be found where the L-curve does not even have an L-shape so that the method 
cannot be used. 



The interactive method 

The speed and versatility of modern digital computers allows us to restore images 
interactively: the user controls the restorations obtained by means of several values 
of the regularisation parameter and, by tuning A, he selects the best restoration on 
the base of his intuition or of the attainment of some specific purpose. 



QCD condensates from r-decay 

data 



Chapter 3 

The theory of r-decays 



The r is the only lepton heavy enough (m,- = 1.777 GeV) to decay not only into 
other leptons, but into final states involving hadrons as well. These decays offer an 
ideal laboratory for the study of strong interactions, including the transition from 
the perturbative to the non-perturbative regime of QCD in the simplest possible 
reaction. This might explain the tremendous efforts ongoing in both theoretical and 
experimental studies of r physics (for a review see |StOO] ). 
The leptonic decays are: 

r ^Ur + W 

\ / = /X or e. (3.1) 

I- + 1^1 

We will write formulas for r = r"; the formulas for r"*" are obtained with obvious 
changes. Besides the decays (13. ip . we have hadronic decays. At the parton level 
these are given by the processes 



Qi + Qj 



(3.2) 



where the fiavour indices i,j stand for the light quarks {u,d,s). The Feynman 
diagram for Eqs. (l3.ip and (13. 2p is shown in Fig J3.1[ 

The permitted hadronic decays may be split into decays that involve only non- 
strange particles 



d + u, 



(3.3) 



or decays involving strange particles 



r ^ Ur + W 



s + u. 



(3.4) 



Although the theory will be presented in general, i.e., for both types of decay, we 
are especially interested in the non-strange ones. 
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3.1 Hadronic r-decays 

The decay of the r lepton into hadrons (see Fig J3.1l) was calculated by Paul Tsai 
|Ts71] . even before the discovery of the r. Ignoring the propagator of the boson, 
the matrix element is given by the product of the leptonic and the hadronic current: 

G 

iM = -^(r I Ji^pJhad,|.| had, z/^), (3.5) 

where is the Fermi constant. The leptonic current is the standard left-handed 
one, Jjgp = f7^(l — 75)z/t-, while the hadronic current can be any of the vector 
•^had = ^ij^iT^^j, axial-vector current J/^^^ = Y.i,j^ij^il''l^<li combinations 
of them. By convention, we have included in the definition of the hadronic current 
part of its coupling to the W boson, the CKM mixing matrix element Vij. Other 
factors were explicitly taken into account when writing out the matrix element (13.51) . 

With this matrix element, one can calculate the decay width T[t Vr + had) 
as [BNP92llNa88a] 



dV{r ^v, + had) = ^(27r)^5W(phad - 0)^1^" 

(3.6) 

x^W rf4xe*<?^(O|Jhad,^(a;)|had)(had|J^,(O)|O)ci(^hadC?0,,, 

had 

where del) denotes the invariant phase space elements of the hadrons and the neutrino, 
and L^'^ is the leptonic tensor. The total 4-momentum of the hadronic system is 
written as q {q^ = s). Now the optical theorem (Fig J3.2p can be used to write the 
matrix element for the production of hadrons in terms of the imaginary part of the 
forward scattering amplitude: 

dTir ^ z/.+had) = ^^L^''^2Im^ J d^x e''^^0\TA^,^{x)jl^^M\0)d<l>..- (3-7) 



3.2. Leptonic r-decays 
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Figure 3.2: The optical theorem relates the production of hadrons from r 
decays (left diagram) to the imaginary part of the forward scattering amplitude 
of the diagram on the right. 



This is quite an important step. While (13. 6p requires the calculation of matrix 
elements of exclusive final states (and their summation), Eq. (l3.7p contains no explicit 
reference to hadronic final states. The first one cannot be handled by perturbative 
QCD, but the second one can. 

We may split the hadron tensor into a transverse and a longitudinal part, writing 

J (3.8) 

= (gMg.-M')n«(?') + ?M9.nW(g2). 

The two functions n*^"'^ introduced here are called the two-point correlators of 
the quark currents (see Section [3^ . They describe the creation of hadronic states 
with total angular momentum J from the vacuum. The integration over the phase 
space of the neutrino can now be performed. The result is 



r(r — > z/t- + had) = 

^ r" il(i_4)^|(i + 24)i,„n"w + i,„n<«>w| 



3.2 Leptonic r-decays 

The leptonic decay width can be calculated in a straightforward way from the Feyn- 
man diagram of Fig J3.1l as well. It is sufficient to treat the process as an effective 
four-fermion interaction and add the effect from the propagator as a correction 
later. The result is 

n2 5 

T{r^uM) = ^{l + ^i). (3.10) 

The quantity A; collects a number of corrections shown in Table 13.11 The table 
also shows corrections which are not present in the Standard Model but will modify 
the decay width if new physics is present. Indeed, a massive neutrino would change 
the phase space of the decays and consequently the decay width. The direct limits 
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Standard Model corrections 



phase space correction due to the finite mass mi of 
the charged daughter lepton |BS88] 



QED radiative corrections [Be58l [PGSQl [PCaM 
[MS881 [RSTTl [5178] 

correction due to the W"^ propagator |LY57] 



271 



mi 
nij. 

25 
T 



TT 



5 \m^r± 



Corrections due to new physics 



neutrino mass |BS88] 

scalar current |St94] 

mixing with a 4th generation v |ST97j 



magnetic dipole moment |Ri97] 



electric dipole moment [Ri97j 



-8 



irij. 
mi 



47] 

mr 
- sin^ e 



K K 

2 ^ To 



10 



Table 3.1: Corrections (contributions to A/) to the leptonic width of the r 
lepton. 



on the mass of Ue and are so strict that these cannot have a measurable impact, 
but the Vr could. Also, mixing with a yet undiscovered, heavy, fourth-generation 
neutrino would reduce the decay width. The authors of |D5T98l I5T97] derived from 
universality limits for: 

m^^ < 38MeV, sin^ ^ < 0.008 (3.11) 

at the 95% CL. 

A charged Higgs boson would introduce a scalar coupling and therefore a non- 
vanishing Michel parameter rj, which in turn modifies the decay width. A limit on 
these bosons has been found |St94^ ISt98] : 

m^i > 2.3tan/3, (3.12) 

at the 95% CL. An estimate of the Michel parameter is |5t94] : 

r/ = 0.01 ±0.05. (3.13) 

An anomalous magnetic moment or an electric dipole moment of the charged 
weak current, usually parametrised with the help of parameters k and k |StOO] . 



3.3. The hadronic branching ratio Rr 
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would modify the decay width. The correction is given in Table 13.11 Limits have 
been derived in |D5T98j . They are (95% CL) 

-0.015 < ft < 0.017, |k| < 0.31. (3.14) 

Nevertheless, the corrections are small (< 4%) and can be safely neglected. 



3.3 The hadronic branching ratio 

Rr is defined as the ratio of the total hadronic decay width r(r — > z/,- + had) to the 
leptonic one: 



Rr 



r(r — i/^- + had) 
r(r ^ Uri^ee) 



(3.15) 



With the earlier computed decay widths of Eqs. (13.91) and (13.101) . the hadronic 
branching ratio becomes: 



Rr = 127T 



ds 



mt 



nit 



mt 



1 + 2^ ) Imn(^)(s) + Imn(°)(s) 



(3.16) 



3.4 Operator Product Expansion (OPE) 

As a consequence of asymptotic freedom the theoretical results obtained from QCD 
can be compared with the experimental situation for so-called hard, i.e. high energy, 
processes: at short distances the effective coupling constant cxg becomes small and 
the interaction can be treated perturbatively. On the other hand, any complete 
theory of the strong interaction must include large-distance dynamics as well. In 
particular quarks interact strongly when forming hadronic bound states. 

A great deal of effort has been made towards the construction of new tools for 
reliable computations in the non-perturbative region of QCD. Most of the efforts 
to obtain quantitative results can be divided into two categories |Pa80] : numerical 
computations and analytic calculations. Numerical computations are usually time 
consuming and require powerful computers (sometimes specially designed). They 
are mainly based on lattice gauge theories ^Ko79t IWi74j and are producing promising 
results IHMPR81] . 

There are several approaches using analytical calculations. For several years 
much effort has been devoted to the search for classical solutions of non-abelian field 
theories jCh78] with the hope that a semiclassical approach may shed some light on 
the underlying quantum world and that classical configurations of fields that make 
the action stationary play an important role in the problem of confinement. 

An interesting new approach based on the operator product expansion was 
opened in 1979 |5VZ79at l5VZ79bt l5VZ79c] . This approach is less fundamental in the 
sense that it does not try to solve the problem of confinement but assumes that con- 
finement exists. In practice the effects of confinement can be described through the 
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use of a few parameters, the so called condensates, and this allows one to investigate 
many hadronic properties. 

As stated above, one of the ingredients of this approach is the operator product 
expansion ^Wi69] . Wilson proposed a short distance operator product expansion of 



the following form 

A{x)BiO) - 5^C„(x)a„(0), (3.17) 

n 

where A and B are local operators. The Cn{x) are C-number functions which can 
have singularities on the light cone of the form [x'^—iex^]~^, p being any real number. 
They can also involve logarithms of x^. In general, the complete expansion involves 
an infinite number of non-singular operators On but to any finite order in x only 
a finite number of these operators contributes. The expansion is valid in the weak 
sense: one must sandwich the product A{x)B{0) between fixed initial and final 
states. Similar expansions exist for time-ordered products or commutators. The 
OPE can be verified explicitly for the free scalar and spinor field theories and for 
renormalised interacting fields to any finite order in perturbation theory. In every 
case they are valid for any elementary or composite local fields. 

The nature of the singularities of the functions Cn{x) is determined, in general, 
by exact and broken symmetries of the theory. The most crucial of these symmetries 
is broken scale invariance. The free scalar and spinor field theories with zero mass 
are exactly scale invariant. Mass terms and renormalizable interactions treated in 
perturbation theory break the symmetry but some remainder of scale invariance still 
governs the behaviour of the singular functions |Wi69] . Exact scale invariance means 
that the field theory is invariant under a one-parameter group of transformations 
f/(A). A local operator (9„(x) transforms as 

U\X)On{x)U{X) = A'*(^")C„(Ax). (3.18) 

In free field theories the constant d{On) is the canonical dimension of the field, i.e. 
[On{x)] = M'^^^"), which can be determined for example from the commutation 
relations. 

In an exactly scale-invariant theory the behaviour of the functions C„(x) is de- 
termined, except for some constants, by scale invariance. Performing a scale trans- 
formation in 

A{x)B{y) - ^a(x-i/)a(2/), (3.19) 



—J. -yM 



we obtain 



X'^i^)+'^(^)A{Xx)B{Xy) - Yl ^«(^^ - Xy)X'^'^"^OM. (3.20) 



" n 



If the fields C„(x) are linearly independent, which can always be arranged, one must 
have 

Cn{Xx - Xy) = A'^(^")-'^(^)-'^(^)C„(x - y). (3.21) 



3.4. Operator Product Expansion (OPE) 
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This equation tells us that — y) must be a homogeneous function of order 

d{On) — d{A) — d{B) in [x — y). This property as well as the behaviour under 
Lorentz transformations determines Cn{x — y) completely, up to a finite number 
of constants |5c71] . When we turn on the interactions |Wi69] . it is still true that 
the C- number functions Cn{x) have a scaling behaviour, as ^ 0, which can be 
summarised with the rules 

Cn{x) x-^", Xn = d{A)+d{B)-d{On), (3.22) 

except that the operator dimensions are no longer given by naive counting of the 
mass dimensions of the operators: they become anomalous in general |Co71] . There 
are, nevertheless, some operators which retain their naive dimensions: these in- 
clude the identity /, and the operators generating symmetries of the theory, such as 
currents J^{x) or the field operator G^'^ . 

The operators On are conveniently classified according to their spin and dimen- 
sion d{On)- In particular, if A and B are scalars, and gauge invariant, only scalar, 
gauge-invariant operators have to be considered. 

Interesting for us is the case of QCD where the two operators A, B are con- 
structed from quark and gluon fields. In this case, the lowest dimension scalar, 
gauge-invariant operators are: 

• dimension 0: / (the unit operator) 

• dimension 4: : m^g^g^ : (no summation in flavours) 

: G'r(x)G^,(x) : 

• dimension 6: : qa{x)Vqa{x)qj3{x)Vqij{x) : 

: qa{x)T{X')apqp{x)q^{x)T{\'')^^q„{x) : 
: mi(f{x)\°-a^uq^{x)G'^''{x) : 
■.fa,cG^''^''{x)G%{x)Gy{x): 

where the first two operators of dimension 6 have been given only for one flavour, 
r denotes any combination of Dirac matrices and therefore for a given flavour there 
are 16 such quantities (J, 7^^, cr'^'^, 75, 75) that render a scalar operator. Any other 
gauge invariant scalar operator of the same dimension can be reduced to these with 
the help of the equations of motion. Notice that the use of the equations of motion 
is legitimate because expectation values in physical states will be taken and the 
equations of motion are satisfied if restricted to the space of physical states. 

Within standard perturbation theory only the unit operator would have a non- 
zero vacuum expectation value, but non-perturbative effects induce non-vanishing 
vacuum expectation values for operators of higher dimensions as well. These are the 
so-called condensates. Therefore non-perturbative effects of QCD introduce power 
corrections of the type > 1, to the perturbative calculation. 
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3.5 Hadronic vacuum polarisation tensor 

The hadronic vacuum polarisation tensor is defined as (cf. Eq. fl3.8l) ): 

%^ = ^J d'x e^''-(0|TJ(;.(a;)J-.(0)t|0), (3.23) 

where J^- = qi'j'^{'j5)qj are the vector (axial-vector) hadronic currents describing 
electroweak interactions, i,j are fiavour indices and T stands for time ordering. 
Performing a Lorentz decomposition, we can split U'^^ into a transverse and a lon- 
gitudinal part (cf. Eq. fl3.8l) ): 

= {q,q. - 9,.q')Tlf^{q') + g,g.<)(g2), (3.24) 
r(^) 



where the indices J = 0, 1 in 11^^, refer to the total spin carried by the correlator. 

t(0) 



The conservation of the vector current implies 11^°^ = 0. The imaginary parts of 



Iljj give the spectral functions: 



l^^lv{s) = -v'^{s), 
Inin«^(s) = ^a^{s), (3.25) 

Im<)^(.) = ^a^i{sl 

which provide the basis for comparing short distance theory with hadronic data. 

Defined as a time ordered product of local operators, the hadronic polarisation 
tensor can be rewritten using the OPE (see Section [3^ : 

where Od = Cd{Od) is the short hand notation for the QCD non-perturbative con- 
densate iPd) of dimension d and its associated perturbative Wilson coefficient Cd'-, 
s = is the momentum transfer. {Od) are vacuum expectation values of gauge 
invariant scalar operators constructed out of quark and gluon fields and called con- 
densates. A more detailed discussion of them is found in Chapter [51 

The dimension = contribution to (13.261) is entirely given by perturbation 
theory (recall that {Oq) = I, Section [231). For that reason it is useful to separate 
the two contributions in fl3.26p 

nS/'-/^.) = n^]ii'^s) + n^)^^(s). (3.27) 

A convenient way to calculate the perturbative part is to use the Adler function 
which is perturbatively constructed as an expansion in the coupling constant: 

D{s) ^ -s-^^m = 4^ E « = V ("^'^ ^ ngp^l^)) • (3.28) 

n>0 



1 



3.5. Hadronic vacuum polarisation tensor 
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The renormalization group equation for a(s), s = —Q^, reads: 

n>0 

da 



d\nQ^ 



(3.29) 



^ 1 

/i2 



a{/.2) 



In the MS scheme for 3 flavours po = 9/4, pi = 4, P2 = 10.06, p3 = 47.23 
|La93[ IRVL97[ ITVZ80] . This allows one to get the perturbative contribution to 
the polarisation operator explicitly at any order of perturbation theory: 



n{s) - n(/x2) 



47r2 



''(^'^ da 

a(At2) P[a) 



(3.30) 



The second term in Eg. (13.27!) is given by the OPE expansion, but now starting 
with dimension d = 2. 



Aim j 




Figure 3.3: Region of validity of the operator product expansion. 



The Adler function constructed in this way has an unphysical cut from s = —Qq 
to s = 0. It is an obvious indication of the fact that QCD is inapplicable at low 
Also, the expression (I3.26P is not valid for all values of s, the series not converging 
for small |s|, where the effective degrees of freedom are hadrons rather than quarks. 
The exact polarisation tensor is known to be an analytical function of s with a 
cut along the positive real semi-axis, while the OPE series with a finite number of 
operators has a very singular behaviour at s = and a possible cut, on the positive 
real semi-axis, starting at some Sq > 0. Based upon these arguments one may draw 
schematically in Fig j3.3l the region of validity of the series (I3.26p . 
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3.6 Dispersion relations 

A dispersion relation connects the real part of an analytic function to its absorba- 
tive (imaginary) part. The name is due to physicists and it is not found in the 
mathematical literature. Similar problems have been studied by the mathemati- 
cians under names such as Hilbert transforms, the Riemann-Hilbert problem, the 
Wiener-Hopf method and Carleman's method. The physicist's name arises from the 
original applications of the technique by Kramers and Kronig to the study of the 
index of refraction describing the dispersion of light waves. More recently, further 
refinements of the method have been employed extensively in elementary particle 
theory as an alternative to perturbative QCD. 

We will first briefly state the mathematical basis without any indication to 
physics. Later on, an application to the two point functions of the previous sec- 
tion will be considered. 

Let us consider a complex-valued function, F{s) (which will later be identified 
with n(s)), of complex argument, s, and assume that: 

• F{s) is real for real s < sq; 

• F{s) has a branch cut for real s > sq', 

• F{s) is analytic for complex s (except along the cut). 



We fix the sign of the absorbative (imaginary) part of F along the branch cut by 



where e > is infinitesimal. Since F is analytic at each point within the contour 



Im s 




Figure 3.4: Contour C of Eq. (l3.32p in the complex s-plane. 



F{s + ie) = ReF(s) + ilmF{s) 



(3.31) 



3.6. Dispersion relations 
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C (Fig J3.4l) . we may apply Caucliy's theorem to find 



27ci Jq s — 
1 / r^' ^^ F{s + ie)- F{s-ie) ^ C F{s) 

TC 



so 



s — — is 2711 /|s|=A2 s — q^ 

where, in the last step, we have used Schwartz' reflection principle: 

F{s + ie)-F{s-ie) = 2i\mF{s). (3.33) 

Suppose that we only know ImF along the branch cut and wish to evaluate F at 
some point q^ . Then, Eq.f l3.32p is not useful since F also appears on the r.h.s. under 
the integral along the circle. However, if 

lim / ds = 0, (3.34) 

then we obtain the unsubtracted dispersion relation 

ne)^'-f\.^^. (3.35) 

Jso s-q^ - IE 

This means that F can be reconstructed at any point from the knowledge of its 
absorbative part along the cut. In particular, the dispersive (real) part of F may be 
evaluated from 

ReF(g2) = Ip r ds ^""^^'^ . (3.36) 
^ Jso s-q^- IS 

In general. Eg. (13.341) is not satisfied, but instead one has 

lim / ds^^^ =a + hq^ + ... , (3.37) 



and thus 



r.A2 

71" 7so s-q^- IS 



1 C ImFfs) 

ReF(g2) = -P / ds- + a + feg^ + ... . (3.38) 



The coefficients of the polynomial in (13.381) . a, 6, depend on the properties of F. 

For the two point functions of the previous section, n(s), one will then have the 
dispersion relation 

I 1 

Ren(g2) = -P / ds Imn(s) + a + hq^ + ... . (3.39) 

Jso "S - q^ 
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Here, the physical meaning of the coefficients a, b, depends of course on the choice 
of the local operator J{x) in the two-point function. In some cases the coefficients 
in question are fixed by low-energy theorems; e.g. if n(0) is known, we can trade 
the constant a in (13.391) for 11(0): 

Ren(g2) = Ren(O) + -P / ds^ + hq^ + ... . (3.40) 

In general, it is always possible to get rid of the polynomial terms by taking an 
appropriate number of derivatives with respect to q^ . 



Chapter 4 

Hadronic spectral functions 



The r lepton was discovered in 1974 from a handful of e — /i events. In the 30 years 
following the discovery, r physics has matured into a field of high precision, testing 
many aspects of our current understanding of particle physics: 



• The coupling of the r lepton to the charged weak current has been tested at 
the level of a few parts per thousand. Its couplings to the weak neutral current 
have been measured with similar precision. 

• Substantial contributions to precision tests of the electroweak theory at LEP 
and SLC have been derived from r production data. 

• The coupling to the photon has been investigated. No anomahes have been 
found so far. 

• More than 100 different, mainly hadronic, branching ratios have been mea- 
sured, allowing tests of QCD and many model predictions. 

• The spectral functions have been determined. They allow one to study per- 
turbative QCD at low energy scales and provide one of the most precise mea- 
surements of the strong coupling constant ag- 

• The structure of the charged current in r decays has been examined for devi- 
ations from V — A. These investigations are sensitive to many kinds of new 
physics. 

• CP violation in the production or decay of the r lepton is a very interesting 
question which has been analysed in a search for a deeper understanding of 
the origin of CP violation. 

• In the light of neutrino oscillations, lepton-flavour-violating r decays are ex- 
pected to occur at some level. 

• The searches for a finite mass of the r neutrino are approaching 10 MeV, 
thereby closing an interesting astrophysical window. 



We are especially interested in the hadronic spectral functions. They are the 
basic ingredients to the theoretical description of r decays, since they represent the 
probability to produce a given hadronic system from the vacuum, as a function 
of its invariant mass squared s. The spectral functions are dominated by known 
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resonances, but tend to approach the quark-level asymptotic regime at m^. They 
have been determined for individual hadronic modes, however they are summed and 
separated into their inclusive vector and axial- vector components for the non-strange 
part. The Cabbibo-suppressed strange part cannot be separated at present, due to 
the lack of necessary experimental information. 

The non-strange r vector spectral functions can be compared to the correspon- 
ding quantities obtained in e~^e~ annihilation by the virtue of isospin symmetry. 
The precision reached in the experimental data makes it necessary to correct for 
isospin-symmetry breaking. Beyond overall agreement, the detailed comparison un- 
veils discrepancies with results from e+e~ data annihilation. The vector spectral 
functions are necessary ingredients to compute vacuum polarisation integrals, re- 
quired for the evaluation of the running coupling constant or the anomalous 
magnetic moment of the muon a^. The disagreement leads to different results when 
using T or e~^e~ spectral functions. While the e"'"e~-based theoretical value of 
disagrees with the measured one by SAa |MRR07] . possibly indicating contributions 
from physics beyond the SM, the r-based calculation is consistent with experiment. 
Although the use of e'^e~ data is a priori more direct, the present situation must 
evolve, requiring more experimental and theoretical cross-checks. 

The observation that hadronic spectral functions derived from r-decay data can 
provide precision information on perturbative QCD is surprising, considering the 
moderate energy scale involved in these decays. Hence it is useful to recall briefly 
the reasons behind this success: 



The total hadronic decay rate normalised to the leptonic width, i?,-, is obtained 
by integrating the total spectral function weighted by a known kinematic factor 
originating from the V — A leptonic tensor. This integral can be transformed 
to a complex contour integral over a circle at |s| = m^, hence involving only 
large complex s values where perturbative QCD can presumably be applied; 

One factor in the weight function, (1 — s/m^)^, cuts off the integrand in the 
vicinity of the real axis where poles of the current correlator are located; 

The invariant-mass spectra in the hadronic r decay, dominated by resonances, 
are thus related to the quark-level contributions expected at large mass. This 
global quark-hadron duality is expected to work best in more inclusive situa- 
tions. In the case of hadronic r decays, this condition seems to be well met, 
with few resonance contributions in each of the equally important vector and 
axial-vector components; 

The perturbative expansion of the spectral function is known to the third 
order in a, |CKT791 ICGSOl [DS79| IGKLQlj [SS9T| : estimates of the fourth order 



coefficient are also available |BCK03| IKS95j ; 



The perturbative and non-perturbative contributions can be treated systema- 
tically using the OPE, giving corrections with inverse powers d of m^; 



4.1. Overview of experiments 
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• The a priori dominant non-perturbative term from the gluon condensate {d = 
4) is suppressed by the second factor of the weight function, (1 + 2s/m^)~^. 
This accidental fact renders hadronic r decays particularly favourable since 
non-perturbative effects are expected to be small; 

• The next d = 6 term is subject to a partial cancellation between the vector 
V and axial-vector A contributions, due to the (1 — 75) factor of the weak 
current. 

4.1 Overview of experiments 

Since its discovery, the r lepton has been studied with ever-increasing precision at 
every new e~^e~ collider that has gone into operation. The appearance of the events 
changed with increasing energy of the machines and improvements in detectors. The 
samples become more numerous and more and more decay modes become available. 

The history of r physics started with a handful of events at SPEAR (MARK 
I), confirmed by the events from DORIS (PLUTO). The energy of DORIS was 
increased in steps to the T(4s) resonance to produce B mesons, but it also was the 
first machine to provide a large number of r pairs, recorded by ARGUS. Meanwhile 
the next generation of accelerators, PEP and PETRA, were built with centre-of-mass 
energies in the continuum around 30 to 40 GeV and with a number of experiments. 
Also CESR started running at the T(4s) with the CLEO experiment, which today 
has the largest r sample. With TRISTAN the Japanese joined, with a machine 
again in the continuum at 50 to 60 GeV. Despite good luminosity the r production 
rate is low, as the cross section falls like 1/s. In 1989 SLC, and shortly afterwards, 
LEP began running at the boson. BEPC was constructed in Beijing to go back 
to the T production threshold and to precisely remeasure, amongst other things, the 
r mass. Table 14.11 lists the accelerators that have produced r pairs and Table 14.21 
summarises the experiments that have analysed them. 

The set of comprehensive measurements of exclusive hadronic branching ratios 
from ALEPH [ALEPHQ^ IALEPHG3] and the non-strange spectral functions from 
ALEPH |ALEPH98[ lALEPHGSj . CLEO |CLE099j and OPAL |OPAL99j have yielded 
important contributions to the measurements of ^^(m^) and study of perturbative 
QCD at low energies. Recent measurements of a set of semi-exclusive branching 
ratios made by DELPHI are also available [DeG4] . 

Although no spectral function measurements have yet been presented by the B- 
factories, BELLE and BABAR, a number of exclusive branching ratio measurements 
of 3-prong and 5-prong final states from BABAR [5o04] are reported. It is evident 
with this data in hand that work on understanding decay mechanisms and form- 
factors for these higher multiplicity states is now required. 

We have chosen to use the final data from the ALEPH collaboration |ALEPH03] 
because, as compared to those available from OPAL |OPAL99] and CLEO |CLE099j . 
they have the smallest experimental errors. Their quality has increased a lot as 
compared to the earlier ones |ALEPH9^ : first, the higher statistics has allowed the 



64 



4. Hadronic spectral functions 



ALEPH collaboration to double the number of bins, and, secondly, the experimen- 
tal errors decreased due to both higher statistics and a better understanding of 
systematic errors. 



Accelerator 


Laboratory 


Energy in GeV 


Years of operation 


SPEAR 


SLAG, USA 


3 - 8 


73 - 88 


DORIS 


DESY, Germany 


8 - 11 


77- 92 


PETRA 


DESY, Germany 


10 - 47 


78 - 86 


PEP 


SLAG, USA 


29 


80 - 90 


CESR 


Gornell, USA 


9 - 12 


79 - 02 


TRISTAN 


KEK, Japan 


50 - 62 


86 - 95 


SLC 


SLAG, USA 


91 


89 - present 


LEP 


GERN, Europe 


88 - 200 


89 - 00 


BEPC 


Beijing, Ghina 


3 - 4 


91 - present 


PEP-II 


SLAG, USA 


11 


99 - present 


KEK-B 


KEK, Japan 


11 


99 - present 



Table 4.1: Accelerators for r physics. 



4.2 Definitions 

The measurement of the r vector and axial-vector current spectral functions requires 
the selection and identification of r decay modes with fixed isotopic G-parity G = +1 
and G = —1, and hence hadronic channels with an even and odd number of neutral 
or charged pions, respectively. Since hadronic final states of different G-parity differ 
also in their quantum numbers, there is no interference between these two states. 
Hence the total hadronic width separates into Fhad = + ^a- 

The spectral functions are obtained by dividing the normalised invariant mass- 
squared distribution dRry/A/ds for a given hadronic mass a/s by the appropriate 
kinematic factor. They are then normalised to the branching fraction of the massless 
leptonic, i.e. electron, channel = (17.810 ± 0.039)% ^ALEPHOSj . 



4.2. Definitions 
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Years 

Experiment Accelerator of 

operation 


i ypicai 

-'-'cm 

m LtC V 


AT , 

pro- 

duced 


Tnt f 

lilt. JL/ 

in 
ill 

pu 


Typical 
e in % 


AT FPN 


T FP 


oy - 


91 


200 


170 


90 


A A/TV 


TPmTAN 

± iViO 


Sfi QA 
ou - cJ^ 


50 - 62 


4 


150 


40 




DORTS 


89 - Q9 


10.58 
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Table 4.2: Experiments in r physics. The number of r pairs produced, 
N^+T--, is given in units of 1000. 'Int. C is the integrated luminosity recorded 
by the experiment. A typical efficiency for the identification of a r pair is 
quoted as e. Only the running periods relevant to the r results are listed. 
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Sew = 1.0198 ±0.0006 accounts for short distance electroweak radiative corrections 
|M588J and the CKM mixing matrix element has the value = 0.9746 ± 0.0006 
[DEHZ03j . Due to the conservation of the vector current, there is no J = contribu- 
tion to the vector spectral function, while the only contribution to Qq is assumed to 
come from the pion pole, oq is connected via partial conservation of the axial-vector 
current (PCAC) to the pion decay constant f-,, through ao,7r(s) = 27i^f^6{s — m^). 

4.3 The mass spectra 

In this section we will sumarise the experimental analysis techniques used by the 
ALEPH Collaboration as described in Ref. |DHZ06j . 

The measurement of the r spectral functions defined in fl4.1l) requires the determi- 
nation of the invariant mass-squared distributions, obtained from the experimental 
distributions after unfolding the effects of measurement distortions. The unfolding 
procedure used by the ALEPH collaboration follows a method based on the regu- 
larised inversion of the simulated detector response matrix using the SVD technique 
|HK96] . The regularisation function minimises the average curvature of the distri- 
bution. The optimal choice of the regularisation strength is found by means of a 
MC simulation where the true distribution is known and chosen to be close to the 
expected physical distribution. 

To measure exclusive spectral functions, individual unfolding procedures with 
specific detector response matrices and regularisation parameters are applied for 
each T decay channel considered. An iterative procedure is followed to correct the 
MC spectral functions used to subtract the cross-feed among the modes. 

Each spectral function is determined in 140 mass-squared bins of equal width 
(0.025 GeV^). 

All systematic uncertainties concerning the decay classification are contained in 
the covariance matrix of the branching ratios. Therefore only the systematic effects 
affecting the shape of the mass-squared distributions, and not its normalisation, 
need to be examined. 

All effects affecting the decay classification and the calculation of the hadronic 
invariant mass are considered in turn. Comparisons of data and MC distributions 
are made and the corresponding biases are corrected, while the uncertainty in the 
correction is taken as input for the calculation of the systematic uncertainty. In the 
case of the spectral functions, the whole analysis including the unfolding procedure 
is repeated, for each systematic effect. This generates new mass distributions which 
are compared bin-by-bin to the nominal ones, hence providing the full 140 x 140 
covariance matrix of the spectral function for the studied effect. 

The systematic studies include the effects from the photon and tt" energy calibra- 
tion and resolution, the photon detection efficiency, the shape of the identification 
probability distribution, the estimate of the number of fake photons, the proximity 
in the calorimeter of other photon showers and of energy deposition by charged par- 
ticles, and the separation between radiative and vr*^ decay photons for residual single 
photons. 
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In addition, systematic errors introduced by the unfolding procedure are tested 
by comparing known, true distributions to their corresponding unfolded ones. 

Finally, systematic errors due to the limited MC statistics and to uncertainties 
in the branching ratios are added. 

In order to illustrate the importance of these systematic uncertainties, one may 
perform an integration over the spectral functions with some given kernel, cha- 
racteristic of a given physical problem. The integration error is then obtained by 
Gaussian error propagation taking into account the correlations. Using moderately 
s-dependent integration kernels, the integration error is dominated by normalisa- 
tion uncertainties, i.e., the errors on the contributing r branching ratios. However, 
the error on an integration with a strongly s-dependent weighting kernel enhancing 
the low energy parts of the spectral functions is dominated by systematics (mainly 
due to the fake photon rejection and the photon efficiency correction at threshold), 
while the central energy region (0.6 - 1.4 GeV^) is statistically limited. When en- 
hancing the higher part of the spectrum, the integration error is equally dominated 
by uncertainties due to the unfolding process, and by limited data and MC statistics. 

4.4 Inclusive non-strange spectral functions 

The ALEPH collaboration |ALEPH05] have provided the spectral functions for the 
hadronic modes, separated into the vector and axial-vector contributions. 

Vector and axial-vector spectral functions 

The inclusive r vector and axial- vector spectral functions are shown in Fig J4.1[ The 
solid line depicts the massless perturbative QCD prediction. Although the statistical 
power of the data is weak near the kinematic limit, the trend of the spectral functions 
clearly indicates that the asymptotic region |GKL91trSS91] is not reached. The QCD 
prediction lies roughly 25% (20%) lower (higher) than the data at for the vector 
(axial-vector) channel. 

To obtain the vector spectral function the two- and four-pion final states were 
measured exclusively, while the six-pion state was only partly measured. The total 
six-pion branching fraction has been determined using isospin symmetry from the 
two- and four-pion ones. 

In complete analogy to the vector spectral function the inclusive axial-vector 
spectral function is obtained by summing up the exclusive axial- vector spectral func- 
tions (odd number of pions in the final state), with the addition of small unmeasured 
modes taken from the MC simulation. 

V zt A spectral functions 

For the total V + A hadronic spectral function one does not have to distinguish the 
properties of the nonstrange hadronic r decay channels. Hence the mixture of all 
contributing non-strange final states is measured inclusively. 
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Figure 4.1: Inclusive vector and axial-vector spectral functions as measured 
by the ALEPH collaboration [ ALEPH05 ]: (a) Vector spectral function vx{s); 
(b) Axial-vector spectral function ai(s). 
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Figure 4.2: Inclusive vector plus axial-vector and vector minus axial-vector 
spectral functions as measured by the ALEPH collaboration [ ALEPHOBj : (a) 
Vector plus axial-vector spectral function v\{^s)^a\{^s); (b) Vector minus axial- 
vector spectral function v\{s) — ai(s). 



The one, two and three-pion final states dominate and their exclusive measure- 
ments are added with proper accounting for (anti-) correlated errors. The remaining 
contributing topologies are treated inclusively, i.e., without separation of the vector 
and axial-vector decay modes. The effect of the feed-through between r final states 
on the invariant mass spectrum is described by MC simulation and thus corrected 
in the data unfolding. In this procedure the simulated mass distributions are itera- 
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tively corrected using the exclusive vector/axial- vector unfolded mass spectra. The 



V + A spectral function is depicted in Fig, 4.2(a) The improvement in precision 



when comparing to the exclusive sum of the two parts in Fig J4.1l is significant at 
higher mass-squared values. 

One nicely identifies an oscillating behaviour of the spectral function due to 
resonance contributions. It does also roughly reach (within errors) the asymptotic 
limit predicted by perturbative QCD at m^. 

In the case of the V—A spectral function, uncertainties on the V/A separation are 
reinforced due to their ant i- correlation. In addition, anti-correlations given between 
r final states with adjacent numbers of pions increase the errors. The V — A spectral 
function is shown in Fig,4.2(b) The oscillating behaviour of the vi and ai spectral 
functions is emphasised but the asymptotic regime seems not to be reached at s = 
m^. However, the strong oscillation generated by the hadron resonances mostly 
averages out to zero, as predicted by perturbative QCD. 



Chapter 5 

Extraction of QCD condensates 



Nowadays, it is reliably established that the microscopic theory of the strong in- 
teraction is QCD, the gauge theory of interacting quarks and gluons. It is also 
established, that unlike, e.g., QED, the vacuum in QCD has a nontrivial structure: 
due to nonperturbative effects, i.e. not admitting the expansion in the intercation 
constant (even if it is small) there persist non-zero fluctuations of gluonic and quark 
fields in the QCD vacuum. The nontrivial structure of QCD manifests itself in the 
presence of vacuum condensates. It is very important to determine these conden- 
sates, to study if one can indeed obtain a consistent description of the low energy 
hadronic physics, to get more insight into the properties of the QCD vacuum, and 
to confront them to theoretical estimates (from instanton calculus or lattice compu- 
tation) in order to test the validity of the condensates from experimental data, but 
also the accuracy of the determination, i.e. to determine their allowed range. 

Condensates, in particular quark and gluonic ones, were investigated starting 
form the 70-ties. Here, first, one should mention the QCD sum rule approach by 
Shifman, Vainshtein and Zakharov |SVZ79at l5VZ79b] , which emphasised the leading 
role of condensates in the calculation of masses of the low-lying hadronic states. 
Starting with this pioneering work, a lot of papers were published on the extraction 
of condensates from experiment. For that purpose one has to relate error affected 
data in the time-like region to asymptotic QCD in the space-like region. This task 
of analytic continuation constitutes, mathematically, an ill-posed problem. In fact, 
extracting condensates from data is highly sensitive to data errors. Not surprisingly, 
results form different collaborations have not been always consistent. 

In what follows we will first discuss general properties of the condensates and 
briefly review previous extractions. Then, a functional approach for the extraction 
of condensates form r-decay data will be presented. 

5.1 Condensates: general properties and previous 
extractions 

As already stated in Section [3751 QCD condensates are vacuum expectation values of 
gauge invariant scalar operators constructed out of quark and gluon fields. We will 
review here the familiar ones, describe their properties and give numerical values 
extracted previously from low-energy and lattice data. 
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5. Extraction of QCD condensates 



Quark condensate: {qq) 

The quark condensate has the lowest dimension {d = 3) and it is one of the 
condensates violating the chiral symmetry. We may rewrite it in the form 

{q<l) = {QLqR + gnqL) , (5.1) 

where qi, qn are the fields of left- and right-chiral quarks. As follows form 
(15.11) . a non-zero value of the quark condensate means the transition of left- 
handed quark fields into right-handed ones and would lead to chiral symmetry 
violation in QCD. (If chiral symmetry is not violated spontaneously, then non- 
zero quark masses will induce (qq) ^ 0.) 

The quark condensate is related by the Gell-Mann-Oakes-Renner relation 
|GMQR68] to measurable constants and the quark masses: 

m = (5-2) 

Here rrin, fn are the mass and the decay constant of the 7r"^-meson (m^r = 
140MeV, = 131MeV), m„ and rud are the masses of u and d-quarks. For 
rUu = 4.2MeV, = 7.5MeV, the quark condensate takes the value: 

(qq) = -1.4 X 10-^GeVl (5.3) 

The value (15. 3p has a characteristic hadronic scale. This shows that chiral sym- 
metry, which is satisfied with a good accuracy in the light quark Lagrangian 
{mu,md/M ~ 0.01, M is the hadronic mass scale, M ~ 0.5 — IGeV), is spon- 
taneously violated in hadronic state spectrum. Due to this fact, the quark 
condensate may be considered as an order parameter for spontaneous chiral 
symmetry breaking. 

Gluon condensate: (— ) 

The gluon condensate has dimesion d = 4 and was initially introduced by Shif- 
man, Vainshtein and Zakharov in Ref . |5VZ79a] and estimated phenomenolog- 
ically by different groups in the literature [BBSlal IBBSlbl iBeSTl [Be88l IBLR85[ 
mm [b05l [KPSM ILNT841 [Ni95l [Ni95l [Ni96l iNiM [Ni02l [NiM IRRY85[ 
l5VZ79allYn99] . Its value is found to be definitely non-zero and positive. How- 
ever, results from different groups differ considerably. Some of the results, 
including the value from Ref. |SVZ79a] , do not satisfy the lower bound derived 
by Bell and Bertlmann |BB81a] from moment sum rules. Recent estimates in 
|Na04] , satisfying this bound from light and heavy quark channels, are of order 
(7.1 ±0.9) X lO-^GeV^ We may add the remark that the gluonic condensate 
conserves chirality. 
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Mixed quark-gluon condensate: g{qG°^J^(T^^q) 

The next in dimension {d = 5) is the mixed quark-gluon condensate which has 
the form: 

9{<lG%^''"^q)^ml{qq). (5.4) 

The value of the parameter ttIq has been estimated from baryon sum rules 
[Bi82l [ChSll [Ch82l [loSTal [loSlbj . the B - B* mass splitting |Na88b] . a string 



model |GS04] and lattice calculations, with a fair agreement. Note that the 
mixed quark-gluon condensate is also violating the chiral symmetry. 

Four-quark condensate: {q^iqq^2q) 

The general form of the dimension d = Q four-quark condensate is as follows: 
{qTiqqT2q), (5.5) 



where Fi, r2 are combinations of Dirac matrices. Following [5VZ79a[ l5VZ79b] 



Eq. (15.51) is usually factorized, and in the sum over intermediate states in all 
channels only the vacuum state is taken into account. The accuracy of such an 
approximation is ~ with Nc the number of colours, i.e., ~ 10%. After 

factorisation, Eq. (15.51) reduces to 

(5.6) 



The four-quark condensate was estimated from different channels |C5SS04l, 
ILNT841 INa04] . The corresponding results indicate a violation of the factorisa- 
tion assumption by at least a factor of 2. 

• Triple gluon condensate: {g^ fabcCG^G^) 

The triple gluon condensate {d = 6) was originally estimated by Shifman, 
Vainshtein and Zakharov |5VZ79a] using an instanton liquid model. Its value 
was confirmed later by lattice calculations [GRBlt RjGPBSj . 

• Higher-order condensates 

The dimension d = 8 and higher dimension condensates were estimated in the 
V |Na95] . -f- A and P + S [NaOlj channels using experimental and lattice 
data [DeOT] . 

More recently, condensates of the V — A channel have been estimated from r- 
decay data [ALEPH98llA5507llBDP506irCGM03l[355M[D5M[Tzmi[N^ 
IRL04t Zy04| and from QCD at large |FGR04j . Large violations of the vacuum 



saturation estimates of these condensates have been observed also in these analyses. 
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5.2 The method: a functional approach 

In what follows we will study a functional methocll], introduced in Ref. |C55S04] . 
which allows, in principle, to extract the condensates within rather general assump- 
tions. 

For our purpose, let us consider a set of functions F{s) which are admissible 
as a representation of the true amplitude if they are real analytic functions in the 
complex s-plane cut along the time-like interval [sq, oo). The asymptotic behaviour 
of F{s) is restricted by fixing the number of subtractions in the dispersion relation 
between F{s) and its imaginary part f{s) = lmF{s + iO) along the cut 

1 

F{s) = — / K{s, z)f{z)dz + suhtr actions . (5.7) 

^ J so 

In general, K{s, z) is the usual Cauchy kernel but we prefer this more general no- 
tation since later on we will use also more complex kernels, e.g. derivatives of the 
Cauchy kernel. 

In order to determine F{s) and f{s) we use the following two available sources 
of information: 

• experimental data measured in the time-like interval Fexp = [sq, Smax], sq > 0: 

/exp('5), 

• theoretical model given by perturbative QCD, i.e., 

— the prediction for F{s) in the space-like interval = [s2,Si]: -Fqcd(s) 

- and fgcnis) = ImFQcois + iO)\^^^^^^^^^.^ since QCD is expected to be 
reliable for large energies. 

As a next step in extracting values for the condensates, we split the integral on 
the r.h.s. of the dispersion relation (I5.7P into two parts: one that can be described 
by the experiment and the other one by the theoretical model, i.e., QCD: 



1 r . . , . , 1 



QCD 



S 



K{s, z)fQCD{z)dz = - / K{s, z)f{z)dz . (5.8) 



QCD prediction: Fqcd{.s) experiment 

The goal of the method is to check if there exists a function F{s) which is in 
accord with both the data on Fexp and the model on F^. For doing this, one can use 
an L^-norni§ approach and define two functionals and xM/] compares 

the true amplitude /(s) with the data. Using its covariance matrix as a weight 
function, we can define: 

1 /* S m ax /* S m ax 

Xl[f] = / dx / dxV-^(a;,x')(/(x)-/exp(a:))(/(x')-/cxp(a;')). (5.9) 

\^ cxpl J So J So 



^The functional method underlying the analysis has first been described in [AM89I IACMQOI ICM90) . 
^TheL^.norm is defined as: = (/,/) = J \f{x)\'^dx. 
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As a measure for the agreement of the true function f{s) with the theory, we 
define by comparing the left and right hand sides of fl5.8p 



XL[f] = m u;l{x) [ Fqcd{x) / K{x,x')f{x')dx' ] dx, (5.10) 



where wl is a weight function for the space-like interval, i.e. an a-priori estimate 
of the accuracy of the QCD predictions and identified with l/aj^{s). o"i(s) is a 
continuous, strictly positive function of s G which should encode errors due to 
the truncation of the perturbative series and the OPE. It is expected to decrease as 
|s| — ^ oo and diverge for s — >■ 0. 

In order to find the true function f{s) one can combine the information contained 
in these two functionals by means of Lagrange multipliers and find the unrestricted 
minimum of 

m = xim + f^xUf], (5.11) 

subject to the condition 

xUf] < Xexp = (5.12) 

which will be the criteria to determine the Lagrange multiplier /z. This procedure 
leads to an integral equation for the imaginary part of the true amphtude, f{x; n) 
[CSS$04| : 

Air I r^'<^^^ r 
f{x]l^)=fe^p{x)-\ PITT-/ dyV{y,x) dx'wL{x')K{x' ,y)FQCD{x') 



+A / dx'K2{x,x')f{x'] fi), 

(5.13) 



so 



where X = 1/ fi and 

|"p I rSmax r 

K2{x,x') = -^^ dyV{y,x) dzwL{z)K{z,y)K{z,x'). (5.14) 

The size of ai determines the minimal value of Xl according to Eq. fl5.10p that 
can be reached by this algorithm. Obviously, widening the error corridor (increasing 
(Tl) will lead to values for Xl mm small as desired. In such a case, the information 
obtained from the fit is not conclusive, since any model function f{s) can be made 
consistent with the data if one allows for a wide enough error corridor. On the other 
hand, narrowing the error corridor will increase ximin) signalling a bad fit, i.e., bad 
consistency of theory with data if the model function is not perfectly describing the 
data. However, a nontrivial result of our approach is the fact, that there exists a 
choice for the error corridor that leads to values ximin = ^(1)- We shall assume 
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that the underlying probabihty distribution is Gaussian, and choose Xl mm 
accordingly 0"^ such that the fit result corresponds to a IcrCL. In practice this is 
done by adjusting the value of the parameters in a^ls). The underlying theory for 
minimisation can be found in Appendix [Dl 

The algorithm to determine acceptable values for the condensates is described 
below. The entire code, programmed in FORTRAN 77, can be found in |AI07] . 



Algorithm 1: Algorithm to determine acceptable values for the condensates. 
Data: fc^p{x), fqcDix), Fqcd{.x), K{x,y). 

Result: values for the condensates Oq and confidence regions around them, 
begin 

initialise M ; /* nr. of free parameters = dimension of Oq */ 
initialise /i ; /* Lagrange multiplier */ 

repeat 

reinitialise fi 

solve Eq. fl5T3|l ^ f{x;fi) 

until xUf] = Xexp 

/opt (a;) ^ 

/^opt 

xlio) ^ xllUtKO) 

ximin ^ minimise xi(^) with respect to O 
solve xliO) = xlmm for O 

solve xliO) = xlmm + for O 
O — Oq defines a confidence region around Oq 
end 



Chapter 6 

V — A analysis 



The V — A correlator, ny^^"* = 11^'' — 11^'' — 11^'', is of particular interest, since 
it contains no truly perturbative contribution in the massless quark limit and it is 
entirely described by the OPE expansion. 

To apply the method described in Section 15.21 to the V — A channel, we identify 
F{s) with Ily^^{s). We will use the appropriate combination of spectral functions, 
Vi{s) — ai{s) — ao(s) as the experimental information. The spin-0 axial spectral 
function, ao{s), is basically saturated by the r rrur channel (see Section and 
can thus be taken into account separately. Inserting this into the dispersion relation 
(15. 7p one finds: 

F{s) = - dzK{s,z)f{z) + i^, (6.1) 



with 

K{s,z) = ^—. (6.2) 
z — s 

Note that for the V — A correlator no subtractions are needed. The term f^/s comes 
from the pion pole; /,r = 0.1307 GeV |PDG06] is the pion decay constant. 
So, one has the information: 

/cxp(s) = 7^ [t;i(s) - ai(s)] , (6.3) 



Fc^cuis) = U^Jis) = Y: 7^ (l + c--«^ + Oial)) (6.4) 

d>6 ^ ^ ^ ^ 

which can be used for the extraction of condensates. The QCD prediction, Fqcnis), 
will then be: 

Fqcd{s) = Fqcd{s) --f-- K{s, z)fQCD{z)dz. (6.5) 

The complete expression for 0^~^ is known to involve two operators [CDGMOU 
ICDGM03j . Assuming vacuum dominance or the factorisation approximation which 
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holds, e.g., in the large- iVc limit, the matrix elements can be written: 

ce + In 

O^o~^ = -ivr«.(gg)2 (f + 327ra,(G„^G'"0) • 

However, our analysis does not rely on the factorisation hypothesis since we are 
going to determine 0^~^, but not the condensates separately. In the result for 
0\f^^ |Zy04| , TTj-Q is defined through the 5-dimensional quark-gluon mixed con- 
densate. Starting from the second order, coefficients in perturbation theory de- 
pend on the regularisation scheme implying that the value of the condensates are 
scheme-dependent quantities. The NLO corrections Cq^'-^ for Oq~^ were computed 
in |LSC86] and the constant cq was found equal to 247/12. This calculation was 
based on the Breitenlohner-Maison definition of 7^ in dimensional regularisation. A 
different treatment of 7^ as used in Ref. |AC94j leads to Cg = 89/12. 

Since the V — A correlator is entirely given by OPE terms and has no truly per- 
turbative contribution, we will take as an error estimate the next higher dimension 
contribution in the OPE, depending on the condensates one wants to determine. 

6.1 1-parameter fit: determination of 0^~^ 

We start with a discussion of results obtained by 1-parameter fits of the dimension 
d = 6 condensate. We have chosen the dimension d = 8 contribution to define the 
error channel on the space-like region, i.e. criis) = \^8~^\uia.^/^'^ with |C^8"~^lm 
in the order of 10~^GeV®. 

6.1.1 Oq~^ at leading-order 

We have performed first the fit for the leading-order (LO) dimension d = 6 conden- 
sate, i.e., we have treated 0\~^ in fl6.6l) as a constant and took c^^*^ = 0. The 
corresponding QCD amplitude, and its imaginary part, are then: 

FQCDis) = ^^-^, fQCDis) = 0. (6.7) 

A typical situation resulting from this algorithm is shown in Fig j6.1l Xl 
indeed a minimum and at least in the vicinity of this minimum it has the expected 



quadratic dependence on Cg (Fig 6.1(a)). The regularised function (Fig 6.1(b)) 



follows nicely the data points, except at large s, as it is well illustrated by Fig,6.1(c 
where the difference between experimental data and the regularised function is plot- 
ted. One can show that the information on the condensate is contained in the lower 
part of the spectrum by adding or removing data-points contained in the saturated 
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upper part: tliis would not cliange tlie result (Fig J6.2l left panel). Thus, the discrep- 
ancies between data and the regularised function at large s, mainly due to larger 
errors in the data, play no role for the determination of O^'"^. 
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Figure 6.1: A typical result of the algorithm: (a) xi [Eg. (15.101) ] as a function 
of Ol'^; (b) the regularised function /(s) [Eq.(l5l3jl] compared with data 
[ALEPH0 5]: (c) difference between data and the regularised function. We have 
chosen for the error parameter O^'"^ = 1.5 x lO^^GeV^. 
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Figure 6.2: 0\~^ as a function of the number of data points used in the 
analysis (e.g. = 120 means that the data points used in the analysis are the 
first 120) (left) and of si, the upper limit of the space-like interval (right). 
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We have restricted the integration range in xt within the hmits S2 < s < si < 0. 
We have observed a well-defined plateau for the value of 0\~^ as a function of 
S\ between —1.5 and GeV^ (Fig |6.2[ right panel) and quote the results for s\ = 
-1.0 GeVl 

We have also studied the behaviour of the algorithm with respect to changes of 
the lower limit, S2- In principle, one could choose S2 to be minus infinity but, as one 
can see in Fig |6.3[ we are forced to choose a finite value. The figure on the right 
shows us that the value of |s2| for which we still have agreement between data and 
experiment is in the order of 10^. For larger values, the minimum of xt (Eq.f lS.lOj) ) 
becomes larger then 1 and thus signalling a bad fit (for more details see Appendix 
ID|) . The reason for this inconsistency, for values of |s2| larger than O(IO^), could lie 
in the assumption of the decoupling of heavy quarks. The heavy quarks could play 
an important role at sufficiently high energies. There are also other reasons, like 
numerical instabilities and rounding errors. However, there exists a plateau for the 
value of 0\~^ as a function of S2 between — lOGeV^ and — 300GeV^. In the rest of 
the analysis we have quoted results for S2 = — ISOGeV^. 




10^ 10^ 10° 
-S2[GeV2 

Figure 6.3: Dependence on the lower end of the space-like interval, S2: Oq~^ 
as a function of S2 (left); XL,mm ^ function of S2, on a log-log scale (right). 



The dimension d = 8 condensate Og~^ was used only to define the error channel 
in the space-like region and thus we expect that the resulting central value of Oq~^ 
does not depend strongly on |C^~^|max5 which is indeed refiected in the plot of 
Fig |6.4l However, increasing |0^~'^|max, i-e., opening the error channel, leads to 
larger uncertainties for Oq~^. A value of xlmm corresponding to a IcxCL can 
be fixed by choosing an error corridor described by the contribution from an Og~^ 
condensate of the expected size of ^ 1.3 x lO^^GeV^. The results of the 1-parameter 
fit at LO can be summarised by quoting the value for 0\~^ |ASSG7] : 

dl-^ = -5.9+11 X lO^^GeV*^ for = 0. (6.8) 
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Tlie value of ^ can be translated into values for the 4-quark condensate 
as{qqf: 



a,{qq? = 2.641°;^^ x IQ-^GeV^ 



(6.9) 



as given by Eg. (16.61) . This value can be compared with the lowest order vacuum 
saturation expression 



<^s{qqY 



2.78 X lO^^GeV^ 



where we have used for the quark condensate the value of Eq. (15.31) . 






O 



V-A 



10° 

[GeV*] X 10^ 



(6.10) 



Figure 6.4: Dependence on the error parameter for the LO analysis: Oq ^ 
as a function of OX~^- 



6.1.2 ^ at next-to-leading-order 

In what follows we are interested in studying the influence of NLO corrections to 
0^~^ and check if the results based on the NLO formula are in agreement with 
the ones found in the last section. For this purpose one needs to consider the 
NLO corrections in (16. 6p and treat 0^~^ as a constant like before. With the NLO 
corrections included, the corresponding QCD amplitude and its imaginary part are: 



Fqcd{s) 



fqcnis) 



(6.11) 



The renormalization scale /i^ was conveniently chosen to be \s\. 

A typical result is seen in Fig J6.5[ Here, one should remark that the discrepancies 
between experimental data and the regularised function are still present. One can 
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thus conclude that including NLO corrections does not improve the quality of our 
fit. The central value of Oq~^ is now shifted with respect to the one from the LO 
fit, Eq. fl6.8p . which was to be expected (Fig J6.6l) . All the consistency checks (see 
Figs. 16.21 and l6.3p performed for the LO fit remain valid. 




Figure 6.5: A typical result of the algorithm for the NLO analysis: (a) 
Xl [Eq. (l5.10p ] as a function of Oq~^; (b) the regularised function /(s) 
[Eg. (15. 13 1 ) ] compared with data [ALEPH05] : (c) difference between data and 
the regularised function. We have chosen for the error parameter Og~^ = 
1.5 X lO^^GeV® and = 89/12. 
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Figure 6.6: Dependence on the error parameter for the NLO analysis with 

C6 = 247/12: C^"^ as a function of O, 
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We can now summarise the results of the NLO fit and quote values for Oq ^ for 
the two available values of the NLO coefficients |A5507] : 

O^-^ = -A.9+li X 10-^GeV'^ for de = 89/12, 

(6.12) 

Og^-^ = -3.6+};^ X 10-^GeV^ for = 247/12. 

The results are based on the 4-loop expression for as with Ajfg{Nf = 3) = 
0.326GeV. They are not sensitive to changing Ajjg within the present experimental 
error of ±0.030GeV. Moreover, the fit results based on the two different values of 
the NLO coefficients ce agree within errors and their difference with respect the 
LO result. Eg. (16.81) . is consistent with a shift calculated from the correction term 
choosing a typical value of 0{1.5)GeV for the renormalization scale in as- This is 
to be expected since our method would work for any s-dependent ansatz for Oq~^ 
as well. 

We have found agreement of theory and data at the IcrCL for the dimension 
d = 8 contribution describing the error with the value ~ 1.3 x 10~^GeV^. 

The values of 0\~^ from Eq. (l6.12p translate into values for the 4-quark conden- 
sate agijiqf according to Eq.( ]6.6p : 

= 2.V9tl% x lO'^GeV^ for = 89/12, 

(6.13) 

as{qqf = l-Sl^o^g x 10"^GeV^ for = 247/12. 

The values are consistent with the one from the LO analysis, Eq.( l6.9p . as well as 
with that of the vacuum saturation expression (I6.10p . 



6.2 2-parameter fit: Oq ^ — Og ^ correlation 

An interesting and new result is obtained from a 2-parameter fit of the dimension 
d = Q and d = 8 condensates. Here we do not include the NLO contribution to Oq~^ 
since the corresponding NLO coefficient for 0^~^ is not known. To define the error 
corridor on the space-like region we have used the dimension d = 10 contribution to 
the operator product expansion and set: 

^L{s) = \Oro~\^^J{-sr (6.14) 

with |Oyo"^l in the order of IQ-^GeV^". 

I I max 

We will not repeat all the consistency checks performed in the 1-parameter case. 
We may argue that they will be qualitatively not changed. Keeping the parameters 
at the quoted values (see Section 16.1.11) . one finds a typical result as the one of 
Fig |6.7[ A direct comparison of experimental data with the regularised function f{s) 



(see Fig 6.7(a)) shows a nice agreement over the full range of s with the exception 
of the highest s-bins. Nevertheless, we have learned from the 1-parameter fits that 
this is due mainly to large experimental errors and that it does not play a sensible 
role in determining the condensates. 
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Figure 6.7: A typical result of the 2-parameter fit (LO): (a) the regularised 
function /(s) [Eq. dS.lSft ] compared with data [ALEPH05] : (b) difference be- 
tween data and the regularised function; (c) 1- and 2a confidence regions in 
the (O^""^, C^~'^)-plane; The central values are marked by dashed lines, (d) 
Xl [Eq. dS.lOp ] as a function of Oq'"^ and Og~^. The error parameter was 
set to \OYn-^\ = 4 X 10-3GeV^°. 
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Figure 6.8: 1- and 2a confidence regions in the (O, 
different values of the error parameter 
marked by dashed lines. 
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The central values are 
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The plot in Fig 6.7(d) shows that Xl ^as the expected quadratic dependence on 
the parameters Oq~ and Og~^, at least in the vicinity of the minimum. However, 
a more convenient way to visualise xi are contour plots, i.e., lines of constant values 
of Xl projected on the (C^""^, -plane. Fig, 6.7(c) shows such contour lines 

corresponding to 1- and 2a confidence regions (CR) for the two parameters (see 
Appendix ID. 21 for the way the confidence regions are defined and constructed). 

As expected, the confidence region becomes larger as we open the error corridor, 
i.e. increase but the central values of Oq~^ and Og~^ corresponding to 

the absolute minimum of xi remain the same within errors. This is well illustrated 
in Fig J6.8[ One can observe that the islands which appear for |C^io~'^|jjja^x = 1 x 
10^^ GeV^° become larger and transform into approximative ellipses as we increase 
l^io~'^lmax sufficiently large |C^io~'^|jjja^x ^^^y become parallelogram like. This 

specific change in the shape has no physical meaning, but ruther shows that at large 
values of 0\~^ and 0\~^ the underlying probability distribution function is not 
Gaussian, x| having a quadratic dependence on the two parameters only in a small 
neighbourhood of its minimum. 
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Figure 6.9: 1- and 2a confidence regions in the 0^~^)-plane from 

the LO 2-parameter fit at the laCL agreement of theory and data. The central 
values are marked by dashed lines. 



The nontrivial result is that we find agreement of theory and data at the laCL 
if we choose |C^io~'^|j„ax ~ ^''^ ^ 10~3GeV^°. The result presented in Fig j6. 91 shows a 
strong correlation of Oq~^ and 0'^~'^. Both the central values as well as the errors 
from the 2-parameter fit are consistent with those from the LO 1-parameter analysis: 
The la range allowed for 0^~^ for fixed 0\~^ = (which is the assumption 
underlying the 1-parameter fit) agrees with fl6.8p and (16.121) . However, leaving the 
value of Og~^ unconstrained, as is the case here, one finds a much larger range for 
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O^-^. The minimum value of Xl located at the values: 
O^-^ = -6.8^2;0 X 10-=^ GeV^ 

(6.15) 

Og^-^ = 3.212-8 X 10-3 GeVl 

The errors on Oq^^ for fixed 0^~^ are small, but the allowed range for 0^~^ is 
not very restrictive (note the different scales for the two coordinate axes in Fig J6.9l) . 
However, the strong correlation allows one to determine a linear combination of 
O^-^ and O^-^ with a rather small error |A55Q7] : 

O^-^ + 2.22 GeV^ C^"^ = -18.30;°:| x lO'^ GeV^ (6.16) 

This is an important result. If one would know one of the parameters from an 
independent determination, one could specify the value of the other one with well 
defined small errors in the order of 3%. 




0^-^[GeV6] X 103 C^-^[GeV6] x 10^ 



Figure 6.10: 1- and 2a confidence regions in the (C^""^, O^-'^j-plane from 
the NLO 2-parameter fit at the IcrCL agreement of theory and data. The 
central values are marked by dashed lines. 

We have also checked that we obtain consistent results when including the NLO 
correction to 0^~^: The 1- and 2a contours are shifted, essentially without changing 
their form (Fig J6.10"l) . to larger values of Oq~^ exactly as can be inferred from 
the 1-parameter fits: for cq = 0(89/12,247/12) we find the minimum of Xl 
O^-^ X lOVGeV^ = -6.8(-6.6,-5.8). 

The linear combinations of Oq~^ and 0\^^ for the two available NLO coeffi- 
cients are then: 

Ol-^ + 2.69 GeV^ Oe^"^ = -14.40^°^^ x lO^^ GeV^ for = 89/12, 

(6.17) 

O^-^ + 4.07 GeV^ O^-^ = -20.80l°|° x lO'^ GeV^ for de = 247/12. 



They are consistent with the LO result and have even smaller errors. 
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6.3 Review and comparison of results 

There exists a number of previous extractions of QCD condensates in the hterature, 
mainly based on sum rule approaches and at LO. They are listed in Table [6TT1 There 
also exist previous extractions at NLO, based on the same functional approach as 
in this work, listed in Table 16.21 

The results for O^'^ cover values between (-2.27 ±0.51) x lO'^ GeV'' |CGM03 
and (— 8.7± 2.3) x 10~^ GeV^ |Na05] . These values correspond to a scale of about 
200 Me V which is comparable to Aqcd- Although errors given by the respective 
authors are typically in the order of 25%, their central values are only in rough 
agreement. The observed variations of these results represent the ambiguities inhe- 
rent in the QCD sum rule approach. Our result nicely falls into the same range, 
also with an error estimate of the same size. 

For previous results range from (-10.8 ± 6.6) x 10"^ GeV^ |BDP506j to 

(15.6 ±4.0) X lO^^GeV^ [Na05] . A recent conservative estimate [RL04] is O^"^ = 
(—121]^^) X 10^'^GeV^. Again, our result agrees within the estimated precision. 
The same is true for 0]'q'^ for which previous results range from (—17.1 ± 4.4) x 
10"=^ GeV^° [Na05] to (78 ± 24) x lO'^ GeV^° [RL04] . The spread of values found in 
the literature for the ci = 12 condensate , Oj^""^, is even larger: from (— 240± 100) x 
10-=^GeV^2 [B DP506] to (14.7 ± 3.7) x lO-^GeV^" [Na05] as one can read off from 
Table 16.11 These values are also consistent with the values we had to choose for the 
error corridor. 

Moreover, it is also interesting to note the agreement of the correlation between 
O^'"^ and 0^~^ with corresponding resuhs from |Na05l |ZyG4] . In Ref. |Na05j . this 
correlation is extracted from weighted finite energy sum rules but no errors are given 
while in Ref. |Zy04] Borel sum rules were used. In the latter, 1-, 2- and 3a confidence 
regions for the correlations OQ~^-Og~^ and Oq~^-OYq~^ are presented. The la 
allowed ranges for C»6^-^ and Cg^"^ are shifted as compared to ours, but the slope 
agrees well within errors. 

Obviously, our numerical results are, from a practical point of view, not superior 
to approaches based on finite energy sum rules. However, the fact that we find 
agreement within errors is not trivial. Since this approach is based on much more 
general assumptions, the results obtained in our analysis give additional confidence 
to the numerical values obtained with the help of QCD sum rules. 
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Table 6.1: Estimated values of the d < 12 ^ condensates in units of 
lO^^GeV^ at leading order. 
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Table 6.2: Estimated values of the d < 8 ^ condensates in units of 
lO^^GeV^ at next-to-leading order. 



^Note the different normalisation of spectral functions. The values shown are adjusted so that they 
can be compared to those from this work. 



Chapter 7 

V ^ A and V + A analysis 



The final precise data on the non-strange hadronic spectral functions from ALEPH 
[ALEPHOS] allows us to perform also an analysis in the V + A, V and A channels 
separately. 

Since the V, A and V + A correlators are dominated by their perturbative contri- 
butions, there are subtractions needed in the dispersion relation (15. 7p . the starting 
point of our analysis. However, as stated in Section 13. 6[ one can eliminate these 
subtractions by taking an appropriate number of derivatives of the correlator in 
question. For V, A and V + A there is only one subtraction and thus one needs to 
take only the first derivative and write 

F{s) = -- [ dxK{s,x)f{x), (7.1) 

^ J so 

where now 

F{s) = -s-^Uv,A,v+Ais), (7.2) 



K{s,x) = ^^^, (7.3) 
[x — s)^ 

and 

fix) = lmUv,A,v+A{x) (7.4) 

is the same as before. 

According to (13.271) . one has: 

F{s) = Dy^AyUs) + E (-sY/2 (1 + OM) ■ (7-5) 



d>4 



The Adler function, -D(s), is known in the massless-quark limit up to terms of 
order at and reads: 



n>0 ^ ^ 



(7.6) 



The coefficients Kn are the same for both A and V channels. For 3 flavours, in 
MS regularisation, Kq = Ki = 1, K2 = 1.64 |CKT79l [CG8OI fPSTO] . K3 = 6.37 
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[GKLQll [559T] and for K4 there are two estimates K4 = 25 ± 25 p<595] and K4 = 
27± 16 |BCK03j . 

The experimental information we use are the appropriate spectral functions, i.e., 
fi(s) for the V channel, ai(s) + ao(s) for the A channel and +ai(s) +ao(s) for 
the V + A one. As seen in Section ao('S) is basically saturated by the pion pole 
and thus we will take it into account separately, as we did before. 

The QCD prediction on the time-like interval (smax, 00) is given by the pertur- 
bative expansion: 

/«"ci = 5^(lH-^ + 0("2)). (T.7) 

To summarise, the QCD information that we have at hand, i.e., the l.h.s. of Eg. (15.81) . 
is: 

= D-(x) + g (1 + cS°'^) . (7,8) 



F^co(^) = fl^(x) + Y: (1 + 4!^^) + k (7.9) 



d>4 



Obviously, for Fq(^^{x) one has to sum the two equations. 

Assuming vacuum dominance, as in the case of the V — A correlator, the ^ + A 
matrix elements can be written as: 



/-OV+A _ '^s r^a \ „NLO _ ^ 

^4: — '^\^ iJLvl } ^4:,V+A — g) 

(7.10) 

OV+A ^ 1^^^ )2 ^NLO = ^ + i! In 4. 

6 s\HH/ , Gy+A 24 18 fi^ 

The coefficients for the a^-corrections were calculated in |CGS85] and |AC94] . re- 
spectively. Considering these matrix elements together with those from the V — A 
channel. Eg. (16.61) . one can easily find similar representations for O^'^. However, as 
before, our analysis does not rely on such representations since we aim to determine 
O^'"^'^^^ and not the condensates (G^j^G^,^), {qq)'^ separately. 

For the error corridor one can use the last known term of the perturbation series 
possibly combined with the first omitted term in the series over condensates. 



7.1 A analysis 

Let us start with the analysis of the axial-vector channel and perform 1- and 2- 
parameter fits in order to determine acceptable values and ranges for the dimension 
d = 4 and d = 6 condensates, as well as their correlation. 



7.1. A analysis 
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7.1.1 1-parameter fit: determination of O 



By a 1-parameter fit we aim to determine the dimension d = A condensate Of at 
leading order. Even though NLO corrections are calculated (see Eq. fl7.10p ). we will 
not take them into account. We have learned from similar fits performed in the 
V — A channel that they are not changing qualitatively the result. 

For defining an estimate for the error corridor in the space-like region one can 
choose to use the last known term of the perturbation series, or the first omitted 
term in the OPE, or a combination of the two of them: 



1 f asi-x) 

J^3 



47^2 



TT 



-x] 



(7.11) 



47r2 



-x) 



TT 



+ 



Oi 



H 2 



-X] 



A typical result of the algorithm is shown in Fig J7.1[ One can see that the 



regularised function (Fig. 7.1(b)) follow nicely the data points, except at large s 



and, as expected, xi the quadratic dependence on Of. 
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Figure 7.1: A typical result of the algorithm: (a) x| as a function of Of; (b) 
the regularised function f{s) [Eg. (15. 131) ] compared with data [A LEPH05 ]: (c) 
difference between data and the regularised function. For the error corridor 
the first omitted term in the OPE was used with \0^\ = 1 x lO^^GeV^ 

I D I max 



We have performed a series of consistency checks summarised in Figs. 17.21 17.31 
and 17.41 i.e., we have studied the behaviour of the algorithm with respect to various 
parameters: the number of experimental data points used in the analysis, the 
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end-points of the time-like interval Tl, si and S2, as well as the dependence on the 
error parameter Oq. 

Like for the V — A channel, one can show that the information on the condensate 
is contained in the lower part of the spectrum by adding or removing data-points 
contained in the saturated upper part of the spectrum: this would not change the 
result (Fig |7.2l . left panel). However, one can observe that there exists also a region 
which deviates from the saturated one, for N < 125. In this region, the high 
oscillation in data points as well as small experimental errors play an important 
role. A good argument to remove these data points from the spectrum, i.e., not 
to use them in the analysis, is to check wether the consistency between theory and 
data is good enough when including them. This is not the case (see Fig J7.2l , right 
panel) and thus for the rest of the analysis we will use only the first = 120 data 
points. 
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Figure 7.2: Of (left) and xlmm ('''ght) as a function of the first data 
points used in the analysis. 
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Studying the behaviour of the algorithm with respect to changes of S2, we found 
that the best simultaneous description of experimental data and theory is when 
we choose to define the error corridor with the help of the last known term in the 
perturbation series, i.e.. 

This is well illustrated in Figs. 17.3] and I7.4( where the dependence of Of and respec- 
tively Xl min "^2 is plottcd for all three possible error corridors. Based on these 
plots we further use the error corridor defined by Eg. (17. 121) and quote results for 
S2 = -3.5GeV2. 

We have also observed a well-defined plateau for the value of Of as a function 
of si between —1.5 and OGeV^ and quote the results for si = — 0.4GeV^. 



7.1. A analysis 
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as a function of the lower end of the space-like interval T^, 



S2, for the three possible error corridors: (a) error corridor defined by the last 
known term in the perturbation expansion, (b) error corridor defined by the first 
omitted contribution in the OPE, (c) error corridor defined by the combination 
of the two previous cases. For the cases (b) and (c), lO^I =3xlO"^GeV^. 




(a) 



(b) 



•4 -2 



(c) 



Figure 7.4: ximm ^ function of the lower end of the space-like interval Tl, 
S2, for the three possible error corridors: (a) error corridor defined by the last 
known term in the perturbation expansion, (b) error corridor defined by the first 
omitted contribution in the OPE, (c) error corridor defined by the combination 
of the two previous cases. For the cases (b) and (c), lO^I = 3x lO^'^GeV^. 



Now, the results of the 1-parameter fit at LO can be summarised by quoting the 
value for O^: 

= 2.8t|]j X 10-^GeV^ at 66.57%CL. (7.13) 

According to Eq. fl7.10p and keeping in mind the chiral symmetry underlying our 
analysis, i.e., assuming O^""^ = 0, one can translate the value of Of into values for 



96 



7. V, A and V + A analysis 



the gluon condensate 

{^G%G%) = 3.361°]^ X 10-^GeV^ (7.14) 

This result is in good agreement with the first extraction by Shifman, Vainshtein 
and Zakharov in |5VZ79a] . 

7.1.2 2-parameter fit: Of — Of correlation 

A two parameter fit is also possible. We aim to fit the condensates of dimension 
d = 4 and c? = 6 at leading order. Based on the results from the 1-parameter fit 
we will take the last known term of the perturbative series as a sensible estimate 
for the error corridor on the space- like region, Eg. (17. 121) . All the consistency checks 
performed for the 1-parameter fit remain valid and thus we keep the corresponding 
parameters, i.e., Si, S2 and A^, at the quoted values (see Section 17.1.11) . 

The result of the 2-parameter analysis is presented in Fig |7.5[ This result cor- 
responds to a 62.58%CL agreement of experimental data with theory. As in the 
1-parameter direct comparison of experimental data with the regularised 



function f{s) (see Figs. 7.5(a) and 7.5(b)) shows a nice agreement over the full 



range of s with the exception of the highest s-bins. We can then conclude that 
including the contribution from the next higher-order condensate in the OPE does 
not improve the quality of the fit. 



Also Xl (Fig 7.5(d)) has the expected quadratic dependence on the two free 



parameters, and O^. Fig 7.5(c) shows lines of constant values of xi projected 
on the {of, (9^) -plane corresponding to 1-, 2- and 3(jCR around the two parameters 
(for a definition of these confidence regions one can consult Appendix ID. 2p . As seen, 
there exists a strong negative correlation of Of and Of, the central values, i.e., the 
values corresponding to xi.min? being located at: 



Of = 3.8tol X 10-3GeV^ 



Of = -1.0^°:^ X lO-^GeV^ 



(7.15) 



These values are at the same 62.58%CL for the agreement between theory and data. 
The errors presented are the la errors. 

The strong correlation allows one to determine a linear combination of Of and 
of with rather small errors: 

of + 0.65GeV2 ■ Of = 1.60l°:| x 10-^GeV^ (7.16) 

and thus if one of the parameters is known from an independent determination, one 
could specify the value of the other one with well defined errors. 

These results are in agreement with those from the 1-parameter fit. However, 
since here the value of Of was left unconstrained, we have found a larger range for 
Of. 
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(d) 



Figure 7.5: A analysis, 2-parameter fit (LO) at 62.58%CL for the agreement 
between theory and data: (a) the regularised function f(s) [Eg. (15. 13 1)] com- 
pared with data [ALEPHOSj : (b) difference between data and the regularised 
function; (c) 1-, 2- and 3a confidence regions in the (Of , (9^)-plane; The 
central values are marked by dashed lines, (d) xi [Eg. (15. 101] ] as a function of 



Of and Of. 
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7.2 V and V + A analysis 

We will now move to the analysis of the V and V + A channels. We will again 
perform 1- and 2-parameter fits and quote their results. 



7.2.1 1-parameter fits 



Let us first start with the 1-parameter fits and aim to determine the dimension 
d = 4 condensates OY and Oj~^^ at leading order. Like for the A channel, since 
the correlators are dominated by the perturbative regime, there are three possible 
choices for the error corridor on the space- like region, i.e.. 



(Tl(x) 



1 



47r2 



a.A-x 



TT 



V 



(7.17) 



47r2 



-X] 



TX 



-\ 2 



+ 



■ OX 

3 ^ 



— x)^ 



for the V channel. For the V+A channel one needs to add the separate contributions 
coming from the V, Eq. fl7.17p . and A, Eq. fl7.11l) . channels. One should keep in mind 
+ for any dimension d. 



that 




1 2 3 4 5 6 7 

OXiGeY^] X 10^ 
(a) 



0.5 1 1.5 2 2.5 3 

s[GeV2] 

(b) 



0.5 1 1.5 2 2.5 3 
(c) 



Figure 7.6: A typical result for the V analysis: (a) Xl ^ function of O, 



V. 



(b) the regularised function f(s) [Eq. (l5.13f )] compared with data |ALEPH05] : 

(c) difference between data and the regularised function. We have chosen for 
the error parameter = 1 x lO^^GeV^. 
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Figure 7.7: A typical result for the V + A analysis: (a) Xl ^ function 
of O^'^'^; (b) the regularised function f{s) [Eq. dS.lSj )] compared with data 
[ALEPH05 ]: (c) difference between data and the regularised function. We have 
chosen for the error parameter O^"^^ = 1 x lO^^GeV^. 



Typical results of the algorithm, for both channels, are shown in Figs. 17.61 and 
17.71 Again, as we have seen in the analysis of the V — A and A channels, the 
regularised function follows nicely the data points, except at large s, and Xl the 
expected quadratic dependence. 

We have performed all the consistency checks discussed in the analysis of the A 
channels and found no justification to change the number of data points N used in 
the analysis or si, the upper limit of the space-like interval Tl- Also, the dependence 
on the lower limit S2 shows again that the best simultaneous description of theory 
and data corresponds to an error corridor defined by the last known term in the 
perturbative series. Unfortunately this agreement is very poor now for the V and 
V + A channels, the corresponding Xl min being in the order of 0(10) while a value 
in the order of 0(1) is expected. 

The results of these 1-parameter fits are: 

OY = 2.6l°| X lO-^GeV^ at 0.28%CL, 

(7.18) 

0^+"^ = X lO-^GeV^ at 9.17%CL. 

As before, these values can be translated into values for the gluon condensate: 
{^Gl^Gl^) = 3.12j:|-°^ X 10"2GeV^ from the V channel, 

(7.19) 

ilT^^G^^u) = 3.24+}:°^ X lO-^GeV^ from the V + A channel. 

Both results, Eqs. 17.181 and 17.191 agree with the corresponding results from the 
A analysis, Eqs. 17.131 and I7.14| respectively. Even though in the present analysis 
there is a poor agreement between theory and data, the extracted values for the 



100 



7. V, A and V + A analysis 



condensates seem to be conclusive. As expected, the values of and are 
almost equal and their sum agrees nicely with the one extracted from the V + A 
channel. 



7.2.2 2-parameter fits 

We have also performed 2-parameter fits to determine the correlation between the 
dimension d = 4 and d = 6 condensates. For this purpose, we have kept all the 
other parameters at the quoted values (see consistency checks for the 1-parameter 
fits) and found the results presented in Figs. 17.81 and 17. 9[ Even though these results 
are as expected, i.e., the regularised function /(s) follows nicely the data points and 
xi has a quadratic behaviour, the problem of having ximin ^^e order 0(10) still 
persits and thus the corresponding CL, taken from Fig JD.3l in Appendix [Dl is very 
small. 

One can see that also in these two channels one finds a strong correlation between 
the two parameters with the central values: 



OX = 6.11?:? X 10-3GeV^ 
= -S.St'oi X 10-3GeV^ 



Ol+^ = 9.9tl 'o X 10-^GeV^ 



0^+^ = -A.2t\i X 10-3GeV^ 



at < 0.01%CL (xi,^i„ = 20.42), (7.20) 



at 0.08%CL iximin = 7.11). (7.21) 



Due to the negative correlation. Figs. 7.8(c) and 7.9(c) contain a lot of informa- 



tion about a specific combination of the two parameters, which allows the determi- 
nation of one of the parameters, with well defined errors, if the other one is known 
from an independent analysis. These combinations are: 



+ 0.65GeV2 . = 0.66t°:| x IQ-^GeV^ 

(7.22) 

/4 — z,.z,u„o.5l 



+ 0.65GeV2 . Or+^ = 2.20^°-^? x lO^^GeV^ 



Despite the poor agreement of theory with data, present in the analysis of the V 
and V + A channels, it is important to remark that the values and allowed ranges 
for the condensates are in good agreement with those found from the A and V — A 
channels. Note that the slope of the correlations in the V, A and V + A channels is 
identical, and thus allowing one to see that the sum of the V and A channels agrees 
with that found in the V + A case. Also the central values agree within errors. For 
the V — A case, the values for both condensates, of dimension d = 4 and d = 6, 
agree with the values found by actually taking the difference between the results 
from the V and A channels separately. 



7.2. V and V + A analysis 
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Figure 7.8: V analysis, 2-parameter fit (LO) at < 0.01%CL (xi,mm = 20.42) 
for the agreement between theory and data: (a) the regularised function f{s) 
[Eq. (l5.13 f)] compared with data [ALEPH 05]: (b) difference between data and 
the regularised function; (c) 1-, 2- and 3a confidence regions in the {O^ , Oq)- 
plane; The central values are marked by dashed lines, (d) Xl [Eg. (15.101) ] as a 



function of and O^. 
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Figure 7.9: V+A analysis, 2-parameter fit (LO) at 0.08%CL (xi,min = 7.11) 
for the agreement between theory and data: (a) the regularised function 
f{s) [Eg. (15. 13 1)] compared with data [ALEPH05] : (b) difference between data 
and the regularised function; (c) 1-, 2- and 3a confidence regions in the 
(O^"*""^, 0^"'''^)-plane; The central values are marked by dashed lines, (d) 
Xl [Eq.dsinD] as a function of C^^^ and 0^+^. 
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7.3 Review and comparison of results 

Besides a consistency check of the resuhs presented in this work, one should also 
compare them to those from other groups. There exist in the literature some previous 
extractions of QCD condensates in the V and A channels. They are based on 
sum rules approaches and the analysis is performed at leading order. However, the 
normalisation of spectral functions is different from ours and there are also other 
factors absorbed in the definition of the condensates. We have translated the results 
so that they can be compared to ours. They can be found in Table 7.1. 





V channel 


A channel 




or 










Be88] 




(1.8-4.8) 


- (12.9 - 8.4) 






[DS88J 


(0.6 - 2.8) 


- (8.1 - 4.1) 


(1.2 - 2.5) 


(4.1 - 7.1) 




DS07 




(1.3-4.8) 


- (19.0 - 5.9) 


(1.3-4.8) 


- (16.5 - 1.3) 



Table 7.1: Estimated ranges for the dimension d = A and d = 6 condensates 
of the V and A channels in units of 10^'^GeV'^ at leading order. Existing 
results from the literature are presented. Note that the normalisation, for all 
of them, was adjusted so that they can be compared to those from this work. 





V channel 


A channel 


V + A channel 




or 




Of 




or+^ 




1-parameter fit 


o c+0.9 
-^•"-0.9 




9 O+0.6 




4+1.7 
"-•■^-1.9 




2-parameter fit 


6.1^?:? 


q q+0.7 
'-'•'-'-0.6 


3-8lo:9 




q n+2.1 

J.J_2.0 


-4 2+^-3 



Table 7.2: Results of this work for the dimension d = A and d = 6 condensates 
in units of lO^^GeV^ at leading order. 



One can remark that the values found in this work (Table 7.2) are consistent 
with those from the literature. The sign of O^, though, disagrees with the vacuum 
saturation approximation, provided Oq is negative, as we have found here. The 
negative sign of Oq disagrees also with the results from |D588] . 

As a conclusion, we can state that the values and ranges found for the QCD 
condensates are all consistent among themselves and with previous extractions found 
in the literature even though the agreement between theory and data is very poor in 
the case of the V and V + A channels. However, one can hope that when eliminating 
some (or all) of the restrictions imposed when performing the analysis the agreement 
between theory and data will improve. 

When analysing all four channels, we have assumed chiral symmetry, decoupling 
of heavy quarks and no duality violation. If the chiral symmetry is broken, there 
are also lower order terms entering the OPE and also mass terms would be present 
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both in the OPE and the perturbative expansion. As a consequence, there will be 
also a perturbative contribution to the V — A correlator. These terms could change 
the results, however remains to be seen if the effect is negligible or not. 

Also the inclusion of heavy quarks could change the results. Heavy quarks are 
expected to play an important role at high energies, i.e., they could play a sensible 
role for the lower end of the space-like interval, S2, as well as for Fexp : (smax, oo) 
where several resonances are known. For a detailed discussion of these resonances 
one can look for example in |EJKV99j and the references therein. 

Duality refers to the fact that, if the approximation n(s) — *• IIope{.s) were exact 
and this substitution carried no error, one can say that the experimental spectral 
function Imn(i:) is dual to IIope{s). The term duality violation, consequently, refers 
to any contribution missed by this substitution. 

As stated above, in our analysis we have assumed no duality violation. However 
it is interesting to check if, and how much, the results would change if one would con- 
sider duality violation contributions. Unfortunately, the theory of duality violation 
in QCD is almost non- existing. Ref. jShOO] was the first to point out the importance 
of duality violations. That the OPE may also miss other important contributions 
has recently been suggested in |Za03j . A model of duality violation can be found in 
^CGPOSj . 



The inverse conductivity problem 



Chapter 8 

Electrical impedance tomography 



Magnetic resonance imaging (MRI), functional MRI (fMRI), X-ray tomography 
or computed tomography (CT), magnetic encephalography (MEG) and ultrasound 
imaging are examples of medical imaging methods that produce a picture of a pa- 
tient's inner organs without intruding into the body. Images coming from different 
systems complement each other since each of these methods measures the distribu- 
tion of some specific physical parameter, such as mass density (CT) or sound speed 
(ultrasound imaging). 




Figure 8.1: Electrodes attached on the boundary of an object for current 
injection and voltage measurement. 

Some of the methods are more harmful to the patient than others: X-ray mea- 
surements expose her to dangerous radiation whereas ultra sound imaging has little 
or no health risks. Another difference between the systems is that some of the ma- 
chines are very expensive. For these reasons it is valuable to have as many different 
imaging methods available as possible. 

Electrical impedance tomography (EIT) is a medical imaging system that pro- 
duces a picture of the electric conductivity distribution inside the patient. For the 
EIT measurement an array of electrodes is attached for instance around the chest 
of the patient. Measurements are done by feeding electric currents through the 
electrodes and measuring the corresponding voltages at the electrodes. This mea- 
surement is repeated with several different choices of current patterns. The resulting 
image is an approximate map of the electric conductivity. 

Electrical properties such as the electrical conductivity a and the electric sus- 
ceptivity e, determine the behaviour of materials under the influence of external 
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electric fields. For example, conductive materials have a high electrical conductivity 
and both direct and alternating currents flow easily through them. Dielectric ma- 
terials have a large electric susceptivity and they allow passage of only alternating 
electric currents. 






Figure 8.2: EIT is a nonlinear problem. 



The mathematical task of interpreting EIT measurements is an inverse boundary- 
value problem called the inverse electric conductivity problem. It is a challenging 
question in the theory of partial differential equations since it is nonlinear and ill- 
posed. Designing a practical algorithm for EIT is hard because non-linearity makes 
the mathematical solution theoretically difficult and ill-posedness means that a small 
error in the practical measurements may cause a large error in the reconstructed 
image. 

Let us consider a bounded, simply connected region Q C M", for n > 2 and, at 
frequency a;, let 7 be the complex admittivity function 

'y{x,uj) = a{x) + iuje{x), (8.1) 

where cr(a?) is called the conductivity and e{x) the susceptivity. 

The electrical impedance is the inverse of 7(3;) and it measures the ratio between 
the electric field and the electric current at location x & Q. Electrical impedance 
tomography is the inverse problem of determining the impedance in the interior of 
Q, given simultaneous measurements of direct or alternating electric currents and 
voltages at the boundary dQ. 

Different materials have different electrical properties, as shown in Tables 18.11 
and 18. 2| so a map of (t{x) and €{x), for x ^ Q, can be used to infer the in- 
ternal structure in Q. Due to this fact, EIT is an imaging tool with important 
applications in vastly different fields such as medicine, geophysics, environmental 
sciences and non-destructive testing of materials. Examples of medical applications 
of EIT are the detection of pulmonary emboh |CIN99l IH5BB871 [Ho93] . monitoring 
of apnoea |ATW90j . monitoring of heart function and blood- fiow |CFGTV79| IICQOj 



and breast cancer detection |CIN99] . In geophysics and environmental sciences, 
EIT can be useful for locating underground mineral deposits [Pa84], detection of 
leaks in underground storage tanks jRDBLR96j and for monitoring fiows of injected 
fiuids into the earth, for the purpose of oil extraction or environmental cleaning 
|RDBLOC93] . Finally, in non-destructive testing, EIT can be used for the detec- 
tion of corrosion |5KV96] and of small defects, such as cracks or voids, in metals 
IAB5V951 IAR981 ICFMV98[ IESIC891 fFV89llSV9l] . 



8.1. The mathematical model 
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Rock or fluid 1/a [Vt cm) 

Marine sand, shale 1-10 

Terrestrial sands, claystone 15-50 

Volcanic rocks, basalt 10-200 

Granite 500-2000 

Limestone dolomite, anhydrite 50-5000 

Chloride water from oil flelds 0.16 

Sulfate water from oil flelds 1.2 



Table 8.1: Resistivity of rocks and fluids |Ke88] . 



Tissue 


l/cr(r2cm) 


e(/iFm-^) 


Lung 


950 


0.22 


Muscle 


760 


0.49 


Liver 


685 


0.49 


Heart 


600 


0.88 


Fat 


>1000 


0.18 



Table 8.2: Electrical properties of biological tissue measured at frequency 
lOkHz [BB84l[Si<57] . 

EIT has been studied extensively in the last two decades and substantial progress 
has been made in both the theoretical and applied aspects of the problem. At the 
same time, EIT remains an area of active research and it continues to pose a variety 
of challenging questions for theoreticians, numerical analysts and experimentalists 
alike. 

8.1 The mathematical model 

Electric and magnetic flelds with harmonic time dependence 

(8.2) 

n{x,t) = Re{H{x,uj)e'^^}, 
satisfy Maxwell's equations in the following form: 

V X H{x,uj) = 'y{x,uj)E{x,uj), 

(8.3) 

V X E{x,uj) = —iLU^{x)H{x,uj), 

where n{x) is the magnetic permeability. EIT operates at low frequencies u, in a 
regime with admittivities 7 and length scales L satisfying tu/i | 7 | <^ 1, such 
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that, after a simple scaling analysis |CIN99] . equations (18. 3p are approximated by 

V X H{x,uj) = 'y{x,u)E{x,uj), 

(8.4) 

V X E{x,u) = 0. 

We define the scalar electric potential (f) and the vector-valued, electric current 
density with harmonic time dependence T{x,t) = Re{j(a:;, u;)e*'^*}, as 

E{x,uj) = -V(f>{x,u), V X H{x,u) = j{x,uj), (8.5) 

such that the first equation in (18. 4p becomes Ohm's law 

j{x) = -'y{x,uj)V(f){x,uj). (8.6) 

Note that the density of dissipated energy, averaged over a period of oscillations. 



— J X(a3,r) ■£{x,T)dT 



= ^[Re{j{x,uj)} ■ Re{E{x,uj)} + lm{j{x,uj)} ■ lm{E{x,u)} 

= ^aix)\V<Pix,uj)\', (8.7) 

must be strictly positive, so we require that 

cx{x) = Re{7(a;, uj)}>m>0, (8.8) 

with m a positive constant. 

If there are no electrical sources within Q, by definition, j is divergence free so 
Ohm's law (18.61) gives the partial differential equation 

V ■ [7(3;, u)V(j){x, uj)] = in fi, (8.9) 

which one should consider together with either Dirichlet boundary conditions 

(f){x, uj) = V{x, uj), for X e dQ, (8.10) 

or Neumann boundary conditions 

dd)(x uj^ 

'y(x,uj)'V(j)(x,uj) ■ n(x) = 'y(x,uj) — = I(x,uj) ai dQ, (S-H) 

on 

such that 

/ I{x,uj)ds{x) = 0, (8.12) 
Jan 

where n{x) is the outer-pointing normal at x G dQ. 

The boundary value problems (I8.9H8.101) or (I8.9H8.11T) for a known function 
'-f{x,uj) in Q and data I{x,uj) or V{x,uj), given for all x G dQ, are referred to 
as continuum, forward mathematical problems for EIT0. 



^Analytic and semi-analytic solutions to a number of forward problems related to the image recon- 
struction problem of EIT in two and three dimensions can be found in [PKL95j . 



8.2. Modelling the electrodes 
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8.2 Modelling the electrodes 

In practice, we do not know the boundary current I{x,uj) for all x G dQ. What 
we actually know are currents flowing along wires which are attached to discrete 
electrodes, which in turn are attached to the boundary dQ |5CI92] . Then, the 
question is how to model the electrodes? 

The gap model approximates the current density by a constant at the surface 
of each electrode and by zero in the gaps between the electrodes. This model is 
appealing because of its simplicity but it is not accurate |IC91j . 

A better choice is the complete model proposed in |SCI92] . Suppose that Ii{uj) 
is the electric current sent through the wire attached to the l-th electrode. At the 
surface Si of this electrode, the normal current density satisfies 

7(x,^)^^^ds(^) = /,M. (8.13) 

In the gaps between the electrodes, we have 

7(x.c.)5^ = 0. (8.14) 

At the contact of Si with dQ, there is in general (often non-negligible) an electro- 
chemical effect which gives rise to a thin, highly resistive layer. This is taken into 
account by the surface impedance zi{x,uj) and 

(f){x, uj) + zi{x, uj)-f{x, u) ^'^'f ' = Vi{lo) for X e Si, 1 = 1,.. .N, (8.15) 

on 

where Vi{uj) is the measured voltage at the l-th electrode. Finally, due to conserva- 
tion of charge one has 

N 

J]JK^) = 0, (8.16) 
1=1 

and by choice of ground, 

N 

J]^(^) = 0. (8.17) 
1=1 

It is proved in |5CI92j that equations (E9]), (l833ll8T71) have a unique solution and 
that they predict experimental data relatively well. However, the inverse problem 
based on this complete model remains essentially unstudied from the theoretical and 
practical points of view. 

In our analysis and the applications we discuss, we will use the approximation 
that the electrodes are point-like and take the current density to be zero in the 
gaps between them. Such a model is not as accurate as the complete one, but good 
enough for our purpose to first find qualitative results. Later on, one can try to 
implement also more complex models, e.g. gap or complete ones. 
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8.3 Formulation of the inverse problem 

In EIT, the admittivity function '-f{x,uj) is unknown and is to be determined from 
simultaneous measurements of boundary voltages V{x,uj) and current densities 
I{x, uj), respectively. There are two maps, the Dirichlet-to- Neumann and Neumann- 
to-Dirichlet maps, which relate V{x,uj) to I{x,uj). These maps depend non-linearly 
on the unknown 'j{x,u!). 

The Dirichlet-to-Neumann map : H^/\dQ) H-^/\dQ) is defined a^ 

A^V(x,uj)=-f(x,uj)^^^^^ foTxeOn, (8.18) 

on 

where V{x, cu) is arbitrary in if^/^ and (f){x, cu) solves the forward problem (18. 9118. 10]) . 
In the static case u; = 0, we have 7(2;, a;) = (j{x) and = A^-. 

The mathematical formulation of EIT, as first posed by Calderon [CaSOj . is as 
follows: Find the admittivity function 'y{x,u) in L°°{rij^, with a strictly positive real 
part cr{x), given the Dirichlet-to-Neumann map A^. 

In practice, it is not advisable to work with the Dirichlet-to-Neumann map. In- 
stead, one uses the Neumann-to-Dirichlet map (A^)~^ which is smooth and therefore 
better behaved for noisy measurements. Nevertheless, in principle, both maps con- 
tain the same information and, usually, the Dirichlet-to-Neumann map is used for 
convenience. 

The Neumann-to-Dirichlet map (A^)^^ : I H^^'^{dVl) is defined on the re- 
stricted space of currents 



I 



\l{x,uj) e H-^''^{d^) such that / I{x,uj)ds{x) = q\ (8.19) 



and, for any I{x,uj) G I, 

{A-yY^I{x,uj) = (p{x,uo) at dVt, (8.20) 

where (f){x,uj) is the solution of the Neumann boundary value problem (18.9118.111) . 
For cij = 0, we have (A^)"-*^ = {K„)~^. We note that (A^)~^ is the generalised inverse, 
as defined in Section I^TTl of A^. 



^H^/^{dfl) are special Sobolev spaces. H^^/^{dn) is the dual of H^/^{dn). For more details one 
can consult for example [GT83]. 

'^L°°{0,) denotes the Banach space of bounded functions on O with the norm 

l|w|L~(n) = ess sup|u|. 

n 
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8.4 Imaging with incomplete, noisy data 




Vp measured 







> n 1 




measured 









Figure 8.3: Illustration of experimental setups for gathering partial data about 
the Neumann-to-Dirichlet map. 



In practice, we do not have full knowledge of the maps (A^)~^ or A^. Instead, 
we have a set of N experiments, where we define an excitation pattern Ie{x,uj) e I 
and we measure the resulting voltage Ve{x,u), at discrete locations Xp G dQ of the 
electrodes, along the boundary. 

In Fig j8.3l we illustrate two experimental setups: An excitation current /g is 
injected (extracted) at some electrodes and the resulting voltage is measured at all 
or some of the electrodes. The first setup is typical of medical applications, where 
one has access to all points of the boundary (i.e. the surface of the body or parts of 
it). The second setup is typical of geophysics applications, where measurements can 
be made at the earth's surface or in some boreholes. In any case, we say that we 
have partial knowledge of the Neumann-to-Dirichlet map and, furthermore, the data 
Ve are usually contaminated with noise. The practical EIT problem is the following: 
Find the admittivity function 'j{x,lj), given a partial, noisy knowledge o/(A^)~^. 



8.5 A brief history of the problem 

Research on the inverse conductivity problem has two aspects: a purely mathema- 
tical one providing answers to the questions of uniqueness and reconstruction and a 
practical approach having as its goal algorithms that perform well on physical data. 

The theoretically oriented line of research was started by CalderoiS [Ci80]. He 
stated the inverse conductivity problem for the first time (in dimensions n >2) and 
gave first results for conductivities close to a constant. Mathematically his work 
continued the tradition of Gel'fand, Levitan and Borg [Bo48[ IGL51] . 

Once Calderon asked his question, new results started to appear. Kohn and 
Vogelius showed that for piecewise analytic conductivities the Dirichlet-to-Neumann 
map uniquely determines the conductivity in dimensions 2 or more [n > 2) |KV84] . 
Sylvester and Uhlmann showed in dimensions n > 3 that if dQ C C°° then A^. 



"^He studied first electrical engineering before changing to mathematics. 
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uniquely determines a E 

The first reconstruction method was published by A. Nachman jNa88] : this 
method is valid in dimensions strictly greater than two (n > 3). 

Meanwhile the case n = 2 remained unsolved. Why is the two-dimensional case 
more difficult than the higher- dimensional ones? A naive argument for this is that 
in n dimensions the Schwartz kernel of Ao- depends on 2{n — 1) variables while a 
is a function of n variables; this means that there is no over-determinacy in two 
dimensions and all the data have to be used for reconstruction. 

Since the inverse conductivity problem is very complicated mathematically, we 
can expect the practical EIT imaging task to be difficult. In addition to the non- 
linearity and ill-posedness of the theoretical problem there are further problems: the 
Dirichlet-to- Neumann map is not given; instead its inverse operator has to be ap- 
proximated from noisy current-to-voltage measurements with finite-precision. Fur- 
thermore, only a finite number of current patterns can be applied through a finite 
number of electrodes. This situation places restrictions on the resolution of the ima- 
ges and requires the physical model of the electrode system to be refined. This leads 
to the notion of distinguishability describing the resolution it is possible to obtain 
in principle with a given system. 

With these restrictions, many algorithms have been developed and tested on real 
data. They fall roughly in three categories: 

• Solve the linearised inverse problem. This approach works for conductivi- 
ties close to constant; examples include back-projection and one-step Newton 
schemes like NOSER. 

• Solve the nonlinear problem by finding iteratively a conductivity distribution 
that matches the physical measurements best in the least squares sense. These 
methods require regular isat ion due to the ill-posedness of the inverse problem. 

• Solve the problem directly. The methods of this category are called direct 
methods as opposed to iterative ones; they aim as well to solve the full non- 
linear problem. 

Some of these ideas combined with elaborate electrical engineering make it possible 
to build reasonably working EIT devices. This has been done in several laboratories 
around the world. 



"'C*^ is the space of functions having k continuous derivatives with a norm defined as: 

k 

ll/llcna) = Esup|/(")(x)|. 

71 — 

Thus, are infinitely derivable functions with continuous derivatives. 



Chapter 9 

The forward problem 



The attention given to the direct problem of EIT is due to the fact that many image 
reconstruction algorithms are iterative and make many comparisons between the 
boundary data predicted by an estimate of the internal conductivity and the data 
that are measured. Also, data produced by solving the direct problem for known 
conductivities is very usefull in checking and studying reconstruction algorithms. In 
the following we discuss method for solving the direct problem often used in practice. 

For an isotropic conductivity distribution cr{x) the potential 4>{x) satisfies the 
equation (see Section ISAl) 



V ■[a{x)V(p{x)] = 0, (9.1) 

where cr{x) and (f){x) are defined in a volume Q bounded by a closed surface dfl. 
The forward problem described here can be formulated as: // the conductivity a 
is known in Q, together with either the potential or the current ad(j)/dn on the 
boundary dVt, find the potential (f) everywhere inside Vt. 
Eq. fl9.1l) can be written in the form 

A0(a;) = -Y{x), (9.2) 

where 

Y{x) = Z^M . V0(a;) = Vlna(a;) • V(j){x). (9.3) 

(J{X) 

A simplification occurs if o"(a;) differs only slightly from a constant conductivity 
distribution 

a{x) = ao + S(7{x). (9.4) 
Then a linear approximation can be used in (19.21) 

A^{x) = YM^ . v</,o(a;) (9.5) 

where (poix) is the solution of (19.11) for a constant conductivity ctq, i.e. of Laplace's 
equation. 

An elegant way to solve the differential equation (19. 2p is by changing it into an 
integral equation by means of the appropriate Green's function (how this is actually 
done see Appendix |B]). Using the Neumann Green's function one gets: 

(j){x) = il^ix) + [ Y{x')GNix,x')d''x', (9.6) 
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where ip{x) is a harmonic function constructed from the measured Neumann bound- 
ary data: 

Hx) = {(p)an + / ^(P{x')GN{x,x')d^''x'. (9.7) 
Jan 

9.1 Methods based on integral equations 

As it was shown in [CPSSOOl ICP503] , one may use the change of variable r = ^/a to 
transform Eq. fl9.ip into an integral equation for the function \E' = T(f). Indeed, using 
this change of variable, Eq. fl9.1l) becomes 

A^(a;) = -V{x)'^{x), where V{x) = (9.8) 

t{x} 

and, with the help of Green's theorem (see Appendix [B]), one finds the desired 
integral equation: 

"^{x) = {^)an+(f -^'^{x')GN{x,x')r-'x'+ [ V{x')^{x')GN{x,x')rx' . (9.9) 
Jan on' 

Once this equation is solved, the potential 0(a;) can be found anywhere in the domain 
by a simple division of \E' by r. 
Another way to solve the forward problem consists of applying the operator 
V^(ln cT(a;)) ■ V^, to ^Mj to find 

Y{x) = V^{\na{x))-'^^i){x)+ I K{x,x')Y {x')(r'x' (9.10) 

Jn 

where K{x,x') = Va;(ln(7(a3)) ■ 'VxG^ix, x'). 

This is an integral equation for Y{x), the Laplacian of (f){x). Once Y{x) is 
known, one can compute 0(a;) by means of a quadrature using again formula (19.61) . 

Since ^xGn{x, x') is divergent, K{x, x') is not a compact operator, that is, the 
integral equation is not of Fredholm type (see Appendix[Al). However, by considering 
its first iteration, that is: 

Y{x) = V^(ln(T(a;)) ■ V^i^^x) + ! K{x,x')V x'{\i^(y{x')) ■ V x'i^{x')<rx' 

Jn 

+ [ K2ix,x')Y{x')rx' (9.11) 
Jn 

with 

K2{x,x')= [ Vx{\na{x))-VxGNix,y)Vy{\na{y))-VyGN{y,x')d^y, (9.12) 
Jn 

one finds a second kind integral equation with a Hilbert-Schmidt kernel (see Ap- 
pendix . A proof of this statement can be found in Ref . |CP504] , where a similar 
approach for solving the forward problem has been presented. The integral equa- 
tion (19. lip can now be solved by means of usual numerical procedures for Fredholm 
equations. 



9.2. Finite element method 
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9.2 Finite element method 

Another way of solving the direct problem is to use the Finite Element Method 
(FEM). But one should first transform the partial differential equation into a vari- 
ational problem. This is simply done by taking the weighted integral on both sides 
of (HI]) 

0=1 w{x)V ■[a{x)V(p{x)]d''x (9.13) 

with w{x) any square integrable, once differentiable, continuous function, i.e. w G 
Hl{Q). Then, by partial integration, one finds an integral condition for the potential 
inside fl 

= - /" a{x)Vw{x)-V(f){x)d"x+ I cx{x)w{x)-^(j)(x)d''-^x. (9.14) 
Jn Jan on 

It is to be remarked that here the Neumann boundary conditions appear naturally. 

In this way, the problem of finding the potential for a known conductivity 
is changed: Find the potential which, for any function w G Hl{Q), satisfies the 
condition: 

a(xWw(x) ■ V(j)(x)d''x = / w(x)a(x)-^(l)(x)d''-^x. (9.15) 

Jdn on 

Such a formulation is also called a weak formulation. It is weak in the sense 
that the local condition for the potential given by the differential equation is now 
changed into an integral one which is definitely less demanding. 

In order to find a numerical solution of Eq. fl9.15p . let us rewrite it in the form of 
a variational problem: 

a{<p,w) = l{w), \/w (9.16) 

with the bilinear form 

a{u,v)= I aVuVvd'^x, (9.17) 
Jn 

and the linear functional 

l{w)= I fwd^~'x, f = a^. (9.18) 
Jan On 

Choosing {ipi, ip2, ...} as a base in the space of the functions w and postulating the 
expansion 

(P = Y,ZkV>k (9.19) 



k 



for the potential, we find the following system of linear equations to be solved: 

^Zka{Lpk,'fi) = K^i)- (9-20) 
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9. The forward problem 



Figure 9.1: Tent functions for the one-dimensional case. 



For the computation of (otherwise time consuming) integrals we use linear finite 
element spaces |PH92] . It is practical to do so since the bases of these spaces are 
formed out of linear tent-functions (Fig lQ.ip and thus, the integrals differ from zero 
only for nearest neighbour nodal points or when they are taken over the same index. 



i^k^, rf"x = ^ \k-i\> 1, (9.21) 
Vv?fcV</?i = V|A;-i|>l. (9.22) 



The matrix inversion is safe since the matrix is positive-definite. Indeed, 

z^Az = ^ ZiAikZk 



i,k 



= a{u, u) > p\ |m| p. 

Once the system fl9.20p is solved, the direct problem is also solved and thus we 
have the potential in the whole domain and we can further use it to check our 
reconstruction algorithms (see next Chapter). 

For numerically solving the forward problem we have used MATLAB |KBOO] . 
since the FEM is already implemented, for two-dimensional domains, via a toolbox 
(PDF- Toolbox). The main steps of this algorithm are described below, while the 



9.2. Finite element method 
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entire code can be found in [AI0 7j. 

Algorithm 2: Algorithm for solving the 2-dimensional forward problem by means 
of FEM in MATLAB. 

Data: Geometry of the problem, boundary conditions, conductivity 
distribution, desired number of electrodes. 

Result: Data files containing the potential and current on the electrodes. 

begin 

solve the system of equations (19.201) for Zk 

by means of Eg. (19. 191) calculate the potential everywhere inside Q 

calculate the data on the boundary dQ 

write out data files 
end 



Chapter 10 

Reconstruction algorithms 



It is well known that the inverse conductivity problem is a highly ill-posed, non- 
linear inverse problem and that the images produced are very sensitive to errors 
which can occur in practice. There has been much interest in determining the 
class of conductivity distributions that can be recovered from the boundary data, 
as well as in the development of related reconstruction algorithms. The interest 
in this problem has been generated by both difficult theoretical challenges and by 
the important medical, geophysical and industrial application of it. Much theore- 
tical work has been related to the approach of Calderon concerning the bijection 
between the conductivity inside the region and the Neumann-to-Dirichlet operator, 
which relates the distribution of the injected currents to the boundary values of the 
induced electrical potential, [CaSOl IKV851 INa88[ ISUBBj . The reconstruction proce- 



dures that have been proposed include a wide range of iterative methods based on 
formulating the inverse problem as a nonlinear optimisation problem. These tech- 
niques are quite demanding computationally particularly when addressing the three 
dimensional problem. This difficulty has encouraged the search for reconstruction 
algorithms which reduce the computational demands either by using some a priori 
information e.g. [BHOOl IBHV031 ICFMV98j or by developing non-iterative proce- 
dures. Some of these methods [BHOOl IBHV03j use a factorisation approach while 
others are based on reformulating the inverse problem in terms of integral equations 
|CP5500l[rP503llCP504j . 



10.1 Reconstruction from a single measurement 

We will present here a reconstruction algorithm based on reformulating the problem 
in terms of integral equations. A feature of this method is that many of the calcu- 
lations involve analytical expressions containing the eigenfunctions of the kernel of 
these integral equations, the computational part being restricted to the introduction 
of data, the numerical evaluation of some of the analytic formulae and the solution 
of a final integral equation. 

However, this method requires the knowledge on the boundary dQ of the domain 
fl not only of the electrical potential (p and its normal derivative d(f)/dn, but also of 
the electrical conductivity a. In geophysics the conductivity on the surface can be 
measured by taking small soil samples, while in medical applications this might be 
achieved using closely spaced electrodes. 

The method consists in the approximate determination (in the sense of a ge- 
neralised solution of inverse problems presented in Section I2TT1) of Y{x) (Eq. (19.21) ) 
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10. Reconstruction algorithms 



by a single simultaneous measurement of the potential and its normal derivative 
d(j)/dn on the boundary dVl. 



10.1.1 The method 

One can use Eq. fl9.6p to determine the potential on the boundary 

i^i^)\.^an + I Y{x')Gr,{x,x')\^^,^<rx' (10.1) 
Jo. 



with 



^i^)\.ean = {4>)dn + jf^ -^,^{x') G^(a., xO|,,g^ t/^V'. (10.2) 

Eq.f llO.li) represents a generalised Neumann-to-Dirichlet mapping. If the potential 
0(03) is measured on the boundary together with the current distribution a{x)^^^^, 
one can determine Yi clS 3i (generalised) solution of the integral equation fllO.ll) . 
One can rewrite this equation in the form 

n^)Le3n = [ Yi^') G^{x,x')\^^,^d^x\ (10.3) 
Jn 



where 



F{^)Lean = -{<P)dn+<P{x)\^^sn- f^^^<Pi^') GN{x,x')\^^,^d--'x' (10.4) 



contains the measured input. Typically a given current pattern would be applied to 
the surface dQ and the resulting surface potential measured. The integral in (110. 4p 
can be evaluated analytically if the applied current pattern is one of the eigenfunc- 
tions of the problem. Equivalently this integral can be viewed as the potential (f)o{x) 
on the surface dfl due to a constant conductivity density 

M^)Ledn = f^^ Q^A^') GM{x,x')\^^,^d--^x' (10.5) 

The data F{x)\^^q^ are defined on a (n — l)-dimensional manifold dVL, while the 
unknown functions Y{x) live in the n-dimensional world VL. As is usual in inverse 
problems, we seek only a generalised solution of fllO.Sp . i.e., a function Y{x) which 
under the mapping of (110.30 produces an F{x)\^^q^ that is closest, in a least-squares 
sense, to the data input and has minimum norm. We rewrite (110.31) in the form 

F = KY, (10.6) 

with K an integral operator. Its domain is the whole region Vt and its range the 
boundary dVt. If the operator K possesses a singular value decomposition (see 
Appendix O) 

00 

KY = aj{Y, vj)uj, (10.7) 
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where {vj} is a complete set of basis functions on Q, {uj} a complete set of basis 
functions on dQ and aj are the singular values of K, the least-squares solution of 
mTB]) will be 

y+ = ^-(F,«>,. (10.8) 

Due to ill-posedness, one needs to regularise the solution. A convenient tool 
to carry out this operation is the so called truncated singular value decomposition 
(Section 12.31) . The regularised solution will then be 

J 1 

Yres = y,-iF,Uj)v,. (10.9) 

j=l J 

Once Yj.cg{x) is known, the potential 4>^eg{x) can be obtained by means of the 
integral 



^rcg 



'x) = ij{x)+ I Y,^g{x')GN{x,x')(r'x'. (10.10) 



To determine the regularised conductivity crrog(a;), one can use the method of 
characteristics to solve the first order partial differential equation (see Eq. fl9.3l) ): 

V (In dreg (a;)) ■ V0rcg(a;) - l^rcg(a3) = 0, X eVL (10.11) 

subject to the known boundary values a{x)\xedQ.- 

10.1.2 An example: the unit disc 

For the unit disc, Eq. (110.31) takes the form 

F{e)= [ £x'K{e;x')Y{x'), (10.12) 
Jn 

where 

F{9) = + 0(r, 9) |,^^ - ^ 0(r', 9') |,,^, G^(l, 9; 1, 9')d9', (10.13) 

and 

oo 2 

K{9,x') = J2^imui{9)viJr',9'), (10.14) 

l,m=l j=l 



with 



(^Im = -T—ClmJlUD- (10.15) 
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The determination of the eigenvalue expansion for the kernel of Eq. fll0.12p . can be 
found in Appendix IB.3[ 

Eq. fll0.14p represents the singular value decomposition of the integral operator 
(Appendix O. The singular functions {uj{6)} form an orthonormal basis on the 
unit circle (the range of K) and the singular functions {^/^(r, 6*)} one on the unit 
disc (the domain of K) (see Appendix [B]) . The singular values aim vanish for large 
m, I which shows the ill-posedness of the problem (a proof of this statement can be 
found in Appendix IA.3p . The least-squares solution of (110.121) reads 

^"■(^^ ^) = E E — f / de'Fi9')uii9')) vjjr, 6). (10.16) 

~1 '^Im \JdQ, / 



l,m=l j- 



An easy way to compute the integrals over the function F{6') is to take advantage 
of the eigenfunction expansion of the Green's function, which on the boundary reads 
(see Appendix [B]) : 

oo 2 ^ 

G^{1, 9;1,9') = Y.Y. W^K^O, (10.17) 
1=1 j=i 



and write for the data function 

oo 2 „27r 

F(^^) = -(0)an + 0(r,^^)|.=i-EEy"'(^) / 

1=1 j=l 

The integration in (110.161) is now straightforward: 



d^{r', 9') 



dr' 



u 



If a' 



)d9'. (10.18) 



r'=l 



JJd9'F{9')ui{9')=<Pi-jIi, (10.19) 
where we have introduced the abbreviations: 

f.27r 

d9(p{r,9)\r=iui{9), (10.20) 







2lT 



V = / d9 



90(r, 9) 



dr 



u\{9). (10.21) 



r=l 

Thus, the least-squares solution of (110.121) will be: 



oo 2 . X 

^) = E E — - 1^1 (10-22) 



l,m=l j=l 

Performing the regularisation, the final result takes the form: 

L M 2 



1 / 1 \ 

y,eg(r, ^) = E E E — - 0)- (10.23) 

1=1 m=l j=l '™ ^ ^ 



An algorithm where we have implemented all these steps was written in MATLAB 
and can be found in |AI07] . 
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10.1.3 Numerical results 

To illustrate the performance of such an algorithm, we present here a numerical 
example. We try to reconstruct a conductivity distribution consisting of a high/low 
conductivity region within a uniform background. An analytical form of such a 
conductivity can be chosen as: 

y) = 1 + aexp[6((x - xf + {y - yf)\ (10.24) 

having a maximum, for a > 0, or a minimum, for a < 0, at {x,y) = {x,y). 
In our numerical simulations we have chosen two sets of parameters {a,b,x,y): 
ai with (1,-1500,0,0.4) and da with (-0.5,-1000,0.5,0). The two exact con- 
ductivities are displayed in the figure below. 




(a) o-5"°^'=i(a;) (b) af°'^''\x) 



Figure 10.1: Density-plot of the two model conductivities considered (see 
text). 



First, one needs to solve the forward problem as in Chapter [91 Let us consider 
the currents 

as input and simulate the measured values of the potential on the boundary. With 
the simulated data we can easily compute the coefficients (pj and needed for 
the regularised solution Freg(r, 9). To fix the regularisation parameters, i.e., the two 
summation limits L, M in Eg. (110.231) . the interactive method was used and we quote 
results for the values L = 10 and M = 2. 

Yj-egi^x) can already be used directly to obtain rough information on the conduc- 
tivity. One can observe in Fig jl0.2l how the effect of the high/low conductivity region 
can be observed on the boundary. The angular information is fine but unfortunately 
there is no radial information present. 

It is still interesting to check whether one can get a better picture when inclu- 
ding the next step in the reconstruction algorithm, i.e., computing the regularised 



an 



sm{m9) 
cos{m9) 



;i0.25) 
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potential: 



N 



(10.26) 



k=l 



where {wk, r^, 6*^} is a set of N quadrature weights and points for the unit disc 
[EnSO] , and: 



^{r, 9) = {<P)an - ^ j^'^ ln(l + - 2r cos{9 - 0')) ^ 



d9'. 



(10.27) 



r'=l 



The actual number of quadrature weights and points used in the reconstruction was 
N = 172. 






(a) Vri^) 



(b) Yr^{x) 



Figure 10.2: Reconstructed Yj.eg{x) for input current sin(6') [Eq. (ll0.25|) 
Reconstructions for the two model conductivities from Fig JlO.ll are shown. 




(a) cj^r^x) (b) ^ri^) 



Figure 10.3: Reconstructed (p^cgix) for input current sin(^^) [Eg. (110.251) ]. 
Reconstructions for the two model conductivities from Fig JlO.ll are shown. 



10.2. Reconstruction from more measurements 
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The reconstructed potential is shown in Fig |10.3[ One can conclude that there is 
no improvement, even worse, in the reconstruction. We believe that by going further 
and performing also the last step we will not gain more. 

Apart from the unsatisfactory reconstruction presented above, the treatment of 
the inverse conductivity problem with a single set of measurements on two dimen- 
sional domains is from practical reasons not advisable. In this case, the information 
one can use to reconstruct the conductivity is a single set of potentials and currents 
on the boundary of the domain. For an ideal case this information would be of 
infinite order, i.e., the potentials and currents wold be known at any point of the 
boundary, but in practice one knows only a finite number of values, corresponding 
to the number of electrodes. In our numerical example we have used 32 electrodes 
(to be in accord with practice), so that a good quality of the reconstruction was, in 
the first place, not to be expected. Besides these practical reasons, the reconstruc- 
tion algorithm has also a mathematical drawback. In Eg. (110.31) the radial part is 
integrated out, i.e., we loose any radial information which might be present in K. 
Eg. (110.31) maps the domain to the boundary, that is it describes a mapping from n 
to {n — 1) dimensions. The reconstruction algorithm needs then to make the inverse 
transition, from {n — 1) to n dimensions, which is in principle impossible without 
additional information, like a model. 

Even though the radial information is lost, the angular information present in 
Freg can be further used as a gualitative model for similar reconstructions, where 
this a priori information is taken into account. Such an algorithm can be found for 
example in |CP504] . 

10.2 Reconstruction from more measurements 

The linearisation of the problem allows us to use the information from more than one 
single set of measurements for the reconstruction of the conductivity distribution. 
Using as input more and different current configurations one gets a system of integral 
eguations which can be easily solved by means of the generalised inverse (see Section 
12. ip . Once the regularisation parameter is fixed, the conductivity distribution is 
reconstructed by simple matrix multiplication. One needs to perform the matrix 
inversion only once since the measured data are not contained in the matrix. 

10.2.1 Reconstruction by linearisation 

Reconstruction of 5a 

The starting point for the reconstruction of 5a is the basic eguation of EIT (cf. 



In the linear approximation, i.e. when the conductivity distribution differs only 
slightly from a constant 



Eg.daU)) 
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one will have for the potential 

(f){x) = (j)o{x) + 6(f){x), (10.30) 

with (f)o{x) the solution of fll0.28p for the constant conductivity ctq. Inserting fll0.29p 
and ffTU^ into (M^^ one finds: 

V [aoV6(l){x)] = -V [6a{x)V(j)o{x)] ■ (10.31) 

The Neumann boundary condition of the problem was already considered when 
computing (j)o{x) (see Eq. fll0.5p ). so that one needs to solve (110. 3ip together with 
d{ao6(f)) / dn = 0. Considering ctq = 1, the solution will be: 

5^ix)= [ GNix,x')V[6a{x')VMx')]d''x'. (10.32) 
Jn 

The integrand of Eq. (110. 32p can be rewritten as 

Gn{x, x')V [6a{x')VMx')] = V [Gn{x, x')5ct{x')VMx')] 

(10.33) 

-VGNix, x')5a{x')VMx')- 

Taking Sa to vanish on the boundary, the first term in (110. 33p can be omitted, since 
it appears as a surface term in (110. 32p . The solution will then become: 

5(f){x) = - [ VGw(a;,a;')V0o(a;')5c^(a;')(^V. (10.34) 
Jn 

Multiplying both sides with the boundary eigenf unctions, Ui, of Gn{x,x') and in- 
tegrating on the boundary, one obtains an integral equation for 5a: 

{S<p,u,)9n = - [ V(l)iix)VMx)S(T{x)d''x. (10.35) 
Jn 

Reconstruction of In a 

Alternatively, the first iteration of integral equation (19. 6p can be used, i.e.. 



4>{x)=M^)+ I GN{x,x')V<po{x')V\na{x')<rx'. (10.36) 
Jn 

Here, the ground was chosen such that (0)ao = 0. As before, one can transform the 
integrand and, assuming In a vanishes on the boundary, Eq. (ll0.36p becomes: 

f)^[x) = - [ VGjv(a;,a;')V0o(a;')lnfT(a;')c/V. (10.37) 
Jn 



X = 0[X 



Note that here no linearisation of the conductivity was performed. 

Taking the scalar product of 6(f) with the boundary eigenfunctions of Neumann 
Green's function, Ui, one finds an integral equation for In a similar to (110. 35p 

{S(l),u^)9n = - [ V0'o(a;)V0o(a3)ln(T(a;)c?"x. (10.38) 
Jn 
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Comparison: da versus In a 

The kernels of the two integral equations fllO.351) and (110.381) are the same, so the 
quantities to be reconstructed will contain the same information. Indeed, in the 
linear approximation one has: 

In a{x) = In (1 + 5a{x)) ^ 5a{x). (10.39) 

However, the reconstruction of Incr is preferred since in this way the positivity of 
the reconstructed conductivity distribution is assured. 



Discretisation 

One can perform a discretisation of ( 110. 38p and of the domain Q: 

{5(P, Ui)an = -Y^ V0^(a;fc) V0o(a;fe) In a{xk)5V{xk), (10.40) 

k=l 

where Np is number of discretisation points Xk- For sets of measured data, 5(f)j, 
and for Nt "test functions" Ui, one finds the system of equations: 

5(l)ij = ^Kijk\aak (10.41) 

k 

with the short-hand notations 
5(t)ij = {S4>j,Ui)on, 

Ki,k = -V(Pi{xk)VMxkWixk), (10.42) 
In CTfc = lna(xk). 



The matrix inversion will be performed by means of the generalised inverse (see 
Section [2111) and the truncated singular value decomposition used as a regularisation 
scheme. In this way, very small singular values will be cut out and will not enter 
in the reconstruction. This is a rather obvious choice since the measurement errors 
are reinforced by small singular values. The regularisation parameter A, also called 
cut-off parameter, will be chosen in such a way that the errors are not reinforced 
but still the reconstruction is close to the true one. 

For numerically solving Eq. fll0.38p . a MATLAB code was written. The important 
steps of this code are summarised below. The entire code can be found in ^AI07j . 
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Algorithm 3: Algorithm for reconstructing the conductivity distribution from 

more than one set of measurements. 

Data: Geometry and data files. 

Result: Reconstructed conductivity distribution. 

begin 

compute 5(j) [Eq. (fT0:i2li ; 
compute K [Eq. lITIM ] 
repeat 
choose A 

compute the generalised inverse of ^ 
calculate Incr = K^5(j) 
calculate a = exp(lno") 
until \a — (Tcxactl < £, £ small 
end 



10.2.2 Numerical results 

As in Section no. 1.2[ we have performed tests for our algorithm on a two-dimensional 
domain, the unit disc. For this specific case, the kernel of the integral equations 
ffT035D and ffTO^SD becomes: 



r 



|i| + b1-2 



cos((|i| - \j\)9) , i,j > or i,j < 0, 



K,^{r,e) = <^ sin((H - IjI)^) , t > 0,j < 0, |j| > H, (10.43) 



7r 



sin((|j|-|z|)^), z<0,j>0,\j\>\t\. 
The "test functions" are (see Appendix [B|) : 

«^W = ^I '"lltP-^n (10.44) 
^ ' V? \ cos(H^), z < 0, ^ ' 

and the potential due to the constant conductivity o"o and boundary condition m^: 

(Pl{r,9) = ^-l^u,{e). (10.45) 

For the tests we have used several conductivity distributions similar to (110.241) : 
(Ti{x) = 1 + 6ai{x), 

a2{x) = l + 6ai{x) + Scr2{x), (10.46) 
0-3(0;) = 1 + 6ai{x) + 5a2{x) + S(T3{x), 
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with 



6ai{x) = exp -1500 (x^ + (y - 0.4)^)' 



6o'2{x) = — O.Sexp 



-2500 {{x - 0.5)^ + {y + 0.2fy 



;i0.47) 



5(73 (x) = 2 exp 



-1000 + 0.3)2 + + 0.2)2)^ 



The forward problem has been solved by means of the finite element method 
described in Section 19.21 The simulated data, i.e. the potential on the boundary, 
contain no further errors, being exact up to computer accuracy. A first test was 
performed using the "exact" potential and later, additional errors were added to 
check the stability of the reconstruction. For the potential 0o due to a constant 
conductivity, Eq. fll0.5p . the theoretically calculated one was used. For the recon- 
struction, the first 10 input currents were used, that is m = 1, 2, 10 in Eg. (110.251) . 
and both sinus and cosine were considered. The number of electrodes was set to 32. 




Figure 10.4: In the upper part, two- and three-dimensional plots of the model 
conductivity ai, Eg. (110. 461) . are displayed. The corresponding reconstructed 
conductivity distribution is shown in the lower part. The regularisation param- 
eter was chosen to be A = 10~^. 
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Figure 10.5: In the upper part, two- and three-dimensional plots of the model 
conductivity cr2, Eg. (110. 460 . are displayed. The corresponding reconstructed 
conductivity distribution is shown in the lower part. The regularisation param- 
eter was chosen to be A = 10^^. 




Figure 10.6: In the upper part, two- and three-dimensional plots of the model 
conductivity era. Eg. (110. 461) . are displayed. The corresponding reconstructed 
conductivity distribution is shown in the lower part. The regularisation param- 
eter was chosen to be A = 10"^. 
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Figure 10.7: Reconstructed conductivity distribution ai for several values of 
the regularisation parameter A. 




Figure 10.8: Reconstructed conductivity distribution 0-3 from error affected 
data: 5% errors (upper part) and 10% errors (lower part). 
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Figs. 110.41 110.51 , and 110.61 show the three test conductivity distributions from 
Eq. fll0.46p and their reconstruction. One can see that the positions of the pertur- 
bations, as well as their height, are relatively well reconstructed. However, since 
small singular values are cut-off, i.e., there is some information missing, one cannot 
reconstruct neither of them exactly. The reconstruction is stable for a well-defined 
interval of the regularisation parameter, A G [10~'', 10"''], as depicted in Fig jl0.7[ 
There exists no rigorous mathematical reasoning for defining this interval, but one 
can find it by looking at the changes in the reconstructed conductivity distribution 
when scanning A. It is obvious that for too small or too large values, the recon- 
struction will become very bad. This is the case since for too small values of A the 
errors on the data will be reinforced, while for too large values the actual signal gets 
suppressed. 

To check further the stability of our method, we have also added additional 
randomly distributed errors on the data. This is realistic since in practice the 
measurements are always accompanied by errors. For big errors the algorithm looses 
its stability, i.e. the conductivity distribution will be not so well reconstructed. 
Indeed, one can observe in Fig J10.8[ lower part, that the reconstructed perturbations 
become flatter as compared to the error free data (Fig llO.61) or to the one where the 
errors on the data are smaller (Fi ^lO.81 upper part). Also, the position of these 
perturbations is imprecisely reconstructed. These changes are due to the fact that 
the regularisation parameter needs to be increased, i.e., one has to cut more singular 
values as in the case of error free data, and hence loose more information. 



Chapter 11 

Reconstructions based on real data 



One medical problem for which knowledge of internal electrical properties would 
be useful is the detection of pulmonary emboli, or blood clots in the lungs. The 
development of pulmonary emboli is a common, and often serious, complication of 
surgery. Unfortunately, at present the diagnosis is rather involved, one possibility 
being the inhalation of radioactive gas in order to determine the ventilated lung 
region. This is followed by injection of a radio-opaque dye or a dissolved radioactive 
substance into a vein to make an image of the blood circulation. The image of the 
circulation in the lungs is compared with the image of the ventilated region; areas 
that are ventilated but not perfused by blood indicate the presence of emboli. 

Another way to determine the location of gas and blood within the body would be 
to map the internal electric conductivity and susceptivity. These electrical properties 
are very different for air, tissue, and blood (see Table 15?^ Chapter[S]); moreover, they 
vary on different time scales. Thus a time-varying map of the electrical properties 
should show lung regions that are ventilated but not perfused by blood. 

Furthermore, determining the presence of pulmonary emboli from a map of the 
electrical properties would have a number of advantages over present techniques. It 
would require no exposure to X-rays or radioactive material. It would be done at 
the bedside, with a relatively small and inexpensive electrical system. 

Information about the internal electrical properties of a body could have many 
other medical uses. Such information could potentially be used for the following: 
monitoring for lung problems such as accumulating fluid or a collapsed lung, nonin- 
vasive monitoring of heart function and blood flow, monitoring for internal bleeding, 
screening for breast cancer, studying emptying of the stomach, studying pelvic fluid 
accumulation as a possible cause of pelvic pain, quantifying severity of premen- 
strual syndrome by determining the amount of intracellular vs. extracellular fluid, 
determining the boundary between dead and living tissue, measuring local inter- 
nal temperature increases associated with hyperthermia treatments, and improving 
electrocardiograms and electroencephalograms. 

To map the electric conductivity inside a body, a tomograph was constructed in 
collaboration with Dr. K.H. Georgi and N. Schuster. The electrical measurements 
were used to reconstruct and display approximate pictures of the electric conduc- 
tivity inside the body, based on the reconstruction algorithm described in Section 

MM 
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11. Reconstructions based on real data 



11.1 The tomograph 

The tomograph consists of three parts: a sensing head, an electronic device to apply 
and measure electric potential and current patterns, and a computer for the image 
reconstruction (Fig lll.ll) . It was designed for mammography applications. For a 
description of this tomograph we will follow Ref. |AHO507j . More details can be 
found in |Az06] . 




Figure 11.1: The tomograph. 

The sensing head has a diameter of 10 cm and consists of 16 large electrodes, 
arranged on the outer ring of a disk. Through these electrodes the current is injected 
and both current and voltages can be measured. There is another set of 64 small 
high-impedance electrodes placed in the interior which can be used to measure 
additional voltages, however, these measurements have not yet been used to solve 
the inverse problem. 




Figure 11.2: The sensing had. 



11.2. Measurements in a test tank 
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The data acquisition device consists of five modules. The first module is a micro- 
controller to facilitate the communication between the measuring device and external 
components via a RS232 serial interface. The second module generates sinusoidal 
voltages of frequencies in the range of 5-50KHz which can be used to drive 16 (or 32) 
current injecting electrodes. The amplitudes can be set positive or negative by 32 
DACs to 16-bit accuracy and can be modulated to any desired amplitude pattern. 
The resulting current at each injection electrode passes through a precision resistor 
and a special operational amplifier to facilitate the simultaneous measurement of the 
current. The voltages on the interior electrodes are measured with 16-bit accuracy 
with the help of the third module. The fourth module serves to read the data and to 
measure the signal by a peak detector via 8 multiplexers of 16 channels each. The 
measured value of the peak detector is subsequently digitalised to 16-bit accuracy. 
In the fifth module, finally, the sign of the modulation is defined. A more detailed 
description of these five modules is found in |Az06] . 

There are — 1 independent measurements for A^ current injecting electrodes. 
In most reconstruction algorithms it is assumed that currents are prescribed and 
voltages are measured. Because of the much simpler electronic implementation, 
voltages are applied and currents are measured on the injectors (voltages are also 
measured in the interior) in our device. It is, however, no problem to convert 
to current-driven data by linear combination of the various voltage-driven data. 
The current-driven data is preferred since, as stated in Section 18.31 the Neumann- 
to-Dirichlet map (Eq. fl8.20p ) is better behaved then the Dirichlet-to-Neumann one 
(Eq. (18.181) ) for noisy measurements. For optimal resolution it is of advantage to 
apply trigonometric voltage or current patterns of different frequencies |HGI88] . 



11.2 Measurements in a test tank 

So far we have no clinical data at our disposal. Instead, we have placed the sensing 
head into the bottom of an appropriate container of large lateral dimensions and 
filled with a conducting liquid. The level of the liquid in the tank has been kept 
very low so as to approximate a two-dimensional situation. Various objects have 
been immersed into the liquid. Measurements have been taken with and without 
the immersed bodies. The latter serve as reference measurements where necessary. 

Fig jll.3l shows such a phantom and the resulting reconstruction |AHOS07] . A 
metal object of roughly 12 x 13mm had been immersed into the tank. The right- 
hand reconstruction has been computed from a measured reference potential (po 
corresponding to a tank with no object immersed. The one on the left has been 
computed to simulate a situation where only "absolute" data, i.e., the potentials 
are available. Here S(j) has been approximated by eliminating the frequency of the 
injected current from the Fourier spectrum of 0. Both reconstructions are fairly 
good, although the one with "absolute" data is only quahtatively correct. Note that 
potentials and currents are only known on the 16 planar electrodes that are clearly 
visible in the photo of the phantom. 
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Figure 11.3: The phantom and its reconstruction. The reconstruction on the 
left is from "absolute" data, i.e., no measured reference potential available, 
while the one on the right is the reconstruction from so called "difference" 
data, i.e., the measured reference potential was used. 
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Figure 11.4: Singular values of the matrix ( 110. 42p . displayed by the '+' sym- 
bols, and the corresponding singular components of the l.h.s. of ( 110.411) . dis- 
played by the 'o' symbols, for the set of data used in the above reconstruction. 



11.3. Measurements on a human chest 
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For the regularisation, the truncated SVD was implemented. To choose a value 
for the cut-off parameter A the interactive method can be used since we exactly 
know what we want to reconstruct. However, in medical applications one does not 
have, in general, a picture of the expected conductivity distribution. For this reason 
we have tried to develop a more mathematically rigorous way to fix A. The form of 
the general solution to inverse problems (cf. Eq. fl2.14p ). 



suggests that one should compare the singular values aj with the singular compo- 
nents {g,Uj)y. As long as the singular values are larger than the singular com- 
ponents, the errors contained in g (data errors) will be suppressed and when the 
singular values become smaller, the errors will be reinforced. This behaviour is a hint 
to choose the cut-off parameter in such a way that the singular values and singular 
components balance each other. In practice, this is done by plotting aj and {g, Uj)y 
against j and choose for A the point where the two of them cross each other. For the 
case of the measurements in a test tank considered in this section, this plot is shown 
in Fig J11.4[ Here one can recognise the usual behaviour of the singular values, i.e., 
they decrease, and the behaviour of the singular components which are expected to 
oscillate. The value of A used in the truncation is indicated by the dashed lines. 
However, as seen in the previous Chapter, there exists a well-defined interval of 
allowed values for the cut-off parameter, centred at the value chosen above. 

11.3 Measurements on a human chest 

To make measurements on patients, a belt of electrodes is placed around their chest. 
It is important that the electrodes are equally spaced, since the reconstruction al- 
gorithm is dependent on their exact position. 

For such measurements, there exist difficulties with the contact impedance. If 
currents are applied and voltages are measured on the same electrodes, the voltage 
measurements will be sensitive to the contact impedance which is formed at the elec- 
trode. To overcome this inconvenience, we have separated the electrodes measuring 
currents and potentials. If this technique is employed, no current flows through 
the voltage measurement electrodes, which means that no voltage is dropped across 
these electrodes, leading to more accurate voltage measurements across the sample. 

The same tomograph was used, only now the sensing head was changed into 
a 16 double-electrode array. The double-electrodes consist of an outer gold-plated 
ring for current measurements and a middle, also gold-plated, pin for potential 
measurements (see Fig Jll.51) . 

Once again, the reconstruction algorithm from Section [10.21 was used to produce 
pictures of the conductivity distribution inside the human body. Here, we have 
used "absolute" data, i.e., Scj) (see Eq. fll0.42l) ) was approximated by ehminating the 
frequency of the injected current from the Fourier spectrum of (f). As a regularisation 





140 
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Figure 11.5: The double-electrode. 



method the truncated singular value decomposition was implemented, with A chosen 
like in the previous section. 

One such reconstruction is shown in Fig lll.6[ The corresponding plot of the sin- 
gular values and singular components is shown in Fig lll.7l where the dashed lines 
correspond to the values chosen for A. One can recognise the heart (high conductiv- 
ity) and the lungs (low conductivity). Note that this is just a qualitative result. Re- 
member that no reference measurement is at hand and that the background constant 
conductivity was approximated to 1, which is not the case in medical applications. 
Other potential sources of error are the inaccuracies in the location of electrodes. 
These inaccuracies may be in the non-uniformity in the spacing between the elec- 
trodes and/or in the shape of the boundary. However, these reconstructions may be 
useful even though the human volunteer is neither two-dimensional nor circular. 

There is still a long way to go until actual medical applications, such as monitor- 
ing for lung problems, are possible. The resolution of the reconstruction should be 
improved, e.g. by using more electrodes and/or a more complex modelling of them 
(see Section ing . Also the time needed for a full cycle of measurements needs to be 
improved. For applications like the so-called "perfusion", i.e., monitoring the blood 
flow to the lung tissue, a synchronisation of the measurements with the heart beat 
is needed. 
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Figure 11.7: Singular values of the matrix (110.421) . displayed by the '+' sym- 
bols, and the corresponding singular components of the l.h.s. of (110.411) . dis- 
played by the 'o' symbols, for the above set of data. 



Conclusions and further work 



QCD condensates from r-decay data 

Although QCD has been with us for three decades, the knowledge of the values 
of the various fundamental or effective parameters of the theory (with the possi- 
ble exception of the coupling constant), such as quark masses and condensates, is 
still astonishingly limited. The precise data on r-decay obtained by the ALEPH 
[ALEPHOS] and OPAL |QPAL99] collaborations have offered an opportunity for new 
studies, which range from an extraction of the strange quark mass \Gu03\ to the 
determination of various condensate parameters. The dimension d = 6 condensate 
was the most studied one ;ALEPH98llA5507llBDP506llCGM03llC55504l [D504l [TZOll 
[Ni05l IOPAL991 [RUM [2^04] . 

We have used a functional method which allows the extraction of condensates 
from experiment within rather general assumptions. It is obvious that it is an 
inverse problem: we want to determine constants of the theory which describe the 
hadronic structure, from indirect measurements. The data used for this purpose 
was the r-decay data from the ALEPH collaboration [ALEPHOS] at LEP. One has 
to relate error affected data in the time-like region to asymptotic QCD in the space- 
like region. This task of analytic continuation constitutes, mathematically, an ill- 
posed problem. In fact, extracting condensates from data is highly sensitive to data 
errors. Not surprisingly, results from different collaborations have not been always 
consistent. 

The method consists in the comparison of the time-like experimental data with 
the asymptotic space-like QCD prediction in an L^-norm manner. Such a method 
must be well defined from a mathematical point of view and physically meaningful. 
One must give explicitly a quantitative estimate of the errors including both expe- 
rimental and theoretical ones (truncation of the perturbative and operator product 
expansions). Then, a set of parameters is said to be admissible (or compatible, or 
consistent) if there exists at least one function with the required analytical proper- 
ties which goes through the L^-error corridors. Answering to this question can be 
reduced to a minimisation process over a finite number of variables and to a sequence 
of Fredholm equations of the second kind, which leads to feasible computations. 

We have performed 1- and 2-parameter fits using the data available for the V — A 
channel. The first one was used for the extraction of the 4-quark condensate {d = 6) 
both at leading and next-to- leading order in a^. The results are in agreement, within 
errors, with those from conventional analyses based on finite energy sum rules. An 
important result is obtained from the 2-parameter fit, where the free parameters are 
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the 4-quark (d = 6) and quark-gluon (d = 8) condensates. The strong correlation 
found allows one to determine a linear combination of the two condensates with a 
rather small error, in the order of 3%. 

We have also performed 1- and 2-parameter fits for the V, A and V + A channels 
separately. Here, the 1-parameter fit estimates the gluon condensate {d — 4) while 
from the 2-parameter fit a correlation between the gluon and 4-quark condensates 
is found. Again, on the basis of this correlation, one can determine combinations of 
the two condensates with rather small errors. 

We would like to continue the research on this line and: 

• Perform 3-parameter fits in all channels. 

• Include contributions omitted in the present analysis like heavy quarks, bro- 
ken chiral symmetry and duality violation. Such contributions may become 
important for a better stabilisation of the algorithm, i.e., for a more stable 

and continuous dependence of the algorithm (and results) on parameters like 
the end points of the space-like interval and the number N of data points 
actually used in the analysis. 

• To improve the precision and test the consistency, we would like to merge the 
results, find overall correlations and allowed ranges; 

• It would be also interesting to compare, by means of the algorithm presented 
in this thesis, the r-decay data with the data from e"'"e~-annihilation and, 
possibly, resolve the discrepancy between them. 

The inverse conductivity problem 

We have developed an approach for solving the inverse conductivity problem based 
on reformulating it in terms of integral equations and applied it to a two-dimensional 
domain, the unit disc, using no a priori information. Unfortunately, the information 
gained on the reconstructed conductivity distribution is restricted to its angular 
dependence (no radial information is present). One can only hope that using some 
a priori information could improve the reconstruction. 

Also an algorithm based on linearisation was presented and tested both on "ex- 
act" (data produced by solving the forward problem) and real data. A promising 
result is that measurements on a human chest, taken with the tomograph con- 
structed in collaboration with Dr. K.H. Georgi and N. Schuster, yield fairly good 
reconstructions. These results, even though they are just qualitative, encouraged us 
in believing that such an algorithm is useful in medical applications like monitoring 
for lung problems and thus intend to continue our research in this direction: 

• We would like to improve the resolution of the reconstruction by, e.g., using 
more electrodes and/or a more complex modelling for them (e.g. the complete 
model); 
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• Improve the time needed for a full measurement cycle and for the reconstruc- 
tions; 

• Synchronise the measurements with the heart beat, so that one can monitor 
the blood flow to the lung tissue. 

Both approaches presented are geometrically constrained by requiring the avail- 
ability of Neumann Green's function for the Laplace equation, which essentially 

restricts us to rectangular, circular, spherical, cylindrical or half-plane geometries. 
It is a challenging idea to remove this constraint and develop similar approaches 
using the free Green's function. 



Appendices 



Appendix A 

Linear integral equations 



Integral equations are equations in which an unknown function appears under one 
or more integral signs. They occur naturally in many fields of mathematical physics. 
They also arise as representation formulas for the solutions of differential equations. 
Indeed, a differential equation can be replaced by an integral equation which in- 
corporates its boundary conditions. Then, each solution of the integral equation 
automatically satisfies these boundary conditions. Integral equations also form one 
of the most useful tools in many branches of pure analysis, such as the theories of 
functional analysis and stochastic processes. 

We shall first give a classification of integral equations and afterwards concentrate 
on Fredholm integral equations which are among the most popular ones. We will 
consider some specific cases of Fredholm equations of the second kind and aim to 
find their solution. For the Fredholm equations of the first kind, only the case of 
symmetric kernels is considered. Based on the solution for this specific case, we 
argue that integral equations of the first kind can have properties characteristic to 
ill-posedness. For a more detailed theory one can consult for example |CH53a[ Wy28] . 

The short theory of integral equations presented here was very useful throughout 
the entire work. In Chapter [5] we introduced an integral equation of the second kind 
which needed to be solved for some specific cases (Chapters [6] and [7]). When we 
studied the kernel of this equation in more detail, it turned out that we deal with 
Fredholm integral equations of the second kind. Later on, in Chapter [HI Fredholm 
equations of second kind arose again when solving the direct problem of EIT. As 
concerning Fredholm equations of the first kind and their ill-posedness, they were 
needed for the reconstruction of conductivities inside a specified domain from a 
single set of measurements (see Section 110.11) . 



A.l Types of integral equations 

Fredholm equations represent one of the most important classes of linear integral 
equations. We shall distinguish between Fredholm equations of the first and sec- 
ond kind. Equations of the second kind are the most interesting and important in 
applications. The simplest formulation for such equations is 

f{x)=u{x)-\! K{x,t)u{t)dt. (A.l) 

J a 

Here, the unknown function u{x) is a function of the real variable x which, along with 
t, is defined on the interval [a,b]. The interval [a,b] may be finite or infinite. The 
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functions K{x, t) and f{x) are known and defined almost everywhere on a < x, t < 6 
and a < X < b, respectively. K{x, t) is called the kernel of (lA.ip and A is a parameter. 
For Fredholm equations, K{x^ t) satisfies 

/ / \K{x,t)fdxdt <oo, (A.2) 

J a J a 

and f{x) satisfies 

\f{x)fdx <oo. (A.3) 

If f[x) = (or, more precisely, if f{x) vanishes almost everywhere on then 
Eq. (]A.l[) is said to be homogeneous, otherwise non-homogeneous. 

Fredholm equations of the first kind differ from the second kind in that the 
isolated term, u{x), in (lA.ip is absent. The simplest form for it is: 

f{x) = \l K{x,t)u{t)dt, (A.4) 

J a 

where K{x,t) and /(x) satisfy Eqs. (]A.2I) and (lA.Sp . respectively. 
If the kernel of the integral equation has the form 

Kix,t) = ^, (A.5) 

where H{x, t) is bounded and < a < 1, Eq. flA.ip will be called an integral equation 
with a weak singularity. 

When the kernel under consideration is of the form 

K{x,t) = ^^, (A.6) 
X — t 

where the numerator A{x, t) is a different iable function of x and t and if the divergent 
integral is interpreted in the sense of its principal value, Eq. (]A.ip is called a singular 
integral equation. 

Volterra equations of the first, homogeneous and second kinds are defined pre- 
cisely as above, except that b = x is the variable upper limit of integration. There 
exist also other classes of integral equations, e.g. equations with convolution kernels, 
Wiener-Hopf equations, dual equations, integral transforms and non-linear integral 
equations. However, we shall consider here only the case of Fredholm equations. 



A.2 Equations of the second kind 

We will start with integral equations of the second kind since they appear to be 
the easiest in the discussion. We will study a few specific cases, as concerning the 
properties of the kernel, for which a compact solution can be found and also give a 
solution by means of Neumann series. 



A. 2. Equations of the second kind 
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A. 2.1 Degenerate kernels 

A kernel, which can be written as a finite sum of products of hnearly independent 
functions of x and hnearly independent functions of t, 

p 

K{x,t) = J2ai{x)bi{t), (A.7) 



i=l 



is called a degenerate kernel. With the help of this representation and exchanging 
summation with integration, Eq. flA.ip becomes: 

p f-b 

f{x) =u{x) - X^ai{x) hi{t)u{t)dt. (A.8) 
i=i 

One can show that the solution of this Fredholm integral equation reduces to solving 
a system of hnear equations. Indeed, if we define q as the integral in (lA.Sp 

b 

hi{t)u{t)dt, (A.9) 

multiply both sides of (lA.Sp by bm{x) and integrate over x from a to b, we find: 

p 

fm = Cm- X^amiCi, m = l,2,...,p, (A.IO) 

i=l 

where 

fm = I bmix)fix)dx, (A.ll) 



and 

r-6 



bm{x)ai{x)dx. (A. 12) 

Eq. (]A.12p is a set of p linear equations in ci, C2, Cp which can be solved by usual 
methods. 



A. 2. 2 Symmetric kernels 

The theory of integral equations can be developed in greater detail if the kernel 
K{x,t) is symmetric, i.e., if it satisfies the relation 

K{x,t) = K{t,x). (A.13) 



Theorem A.l Every continuous symmetric kernel that does not vanish identically 
possesses eigenvalues and eigenf unctions; their number is denumerably infinite if 
and only if the kernel is non-degenerate. All the eigenvalues of a real symmetric 
kernel are real. 
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Proof. The proof is made in several steps: 

1. Prove the existence of an eigenvalue of a symmetric kernel; 

2. Consider the totality of eigenfunctions and eigenvalues; 

3. Finally prove that all the eigenvalues are real. 

For further details see |CH53a] . □ 

Theorem A. 2 Every continuous function g{x) which is an integral transform with 
symmetric kernel K{x, t) of a piecewise continuous function h{t) 



g{x) = / K{x,t)h{t)dt, (A.14) 

J a 

can be expanded in a series in the eigenfunctions of K{x,t) 

9{x) = ^gi<^i{x), gi = {g,<^i) = (A.15) 

this series converges uniformly and absolutely. □ 
Setting now hiy) = K{y,t), for the 'iterated kernel' 

ir(2)(x,t)= / K{x,y)K{y,t)dy (A.16) 

J a 

we then have the expansion 

K^^)(x,t) = J2^ K{y,t)^.{y)dy (A.17) 
i=i J- 

or 

ir(^)(a:,t) = V^iM?^ (A. 18) 
In the same way, the subsequent iterated kernels 

ir(")(x,t)= [\^''-'\x,y)K{y,t)dy (A.19) 



a 



admit the expansions 



K(-\x,t) = f^^^^^ (n = 2,3,...), (A.20) 



i=l 



all of which converge absolutely and uniformly both in x and in t, and uniformly in 
X and t together. 



A. 2. Equations of the second kind 
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However, such an expansion for the full kernel, 
Ki,,t)^±^^^^, (A,21) 

1=1 

does not always exist. 

We can now derive the formula for the solution of the inhomogeneous integral 
equation flA.ll) for a symmetric kernel K{x, t) which admits the eigenfunction ex- 
pansion flA.2ip . We assume initially that the parameter A is not equal to any of the 



eigenvalues Aj. If the continuous function u with the expansion coefficients {u, (pi) 
were a solution of the integral equation, then, by the expansion theorem applied 
to Xu{t), the function u{x) — f{x) = g{x) would be given by the uniformly and 
absolutely convergent series 

oo 

g{x) = u{x) - f{x) = '^CiLpi{x) = A / K{x,t)u{t)dt, (A.22) 

i=i 

where 

l-b r-b AAA 

Ci = (g, ^Pi) = \ I I K{x, t)ipi{x)u{t)dxdt = — (v9i, u) = — (v^j, /) + — (v^i, g) 

Ja Ja \ \ 

(A.23) 

from which it follows that 

C^ = f^J^, {f^={^^J)). (A.24) 

We thus obtain the series expansion for u, 

oo „ 

u{x) = fix) + XJ2 (A.25) 

1=1 * 

which must represent the solution of (lA.ip . This solution fails only if A = Aj is an 
eigenvalue; it remains valid even in this case if f{x) is orthogonal to all eigenfunctions 
(fi belonging to the eigenvalue A^. 



A. 2. 3 Neumann series and the reciprocal kernel 

The theorems of integral equations given above imply a prescription for actually 
computing solutions with arbitrary accuracy. It does not give the solutions, however, 
in an elegant closed form as in the previous section. 

We write the integral equation (lA.ip . inserting in place of u{t) in the integral 
the expression obtained from flA.ll) . and iterate this procedure. With the aid of the 
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iterated kernels we thus write (lA.ip in the form 

u{x) =f{x) + X K{x,t)f{t)dt +AM K^'^\x,t)u{t)dt 

J a J a 

pb pb 

= f{x) + X K{x,t)f{t)dt +AM K^^\x,t)f{t)dt 

J a 5 

/ K^^\x,t)u{t)dt 



and see that the solution is given by the infinite series 

pb pb 

u{x)=f{x) + X K{x,t)f{t)dt + \^ I K^^\x,t)fit)dt + (A.26) 

J a J a 

provided that this series converges uniformly. If, moreover, the series 

}C{x,t) = K{x,t) + XK^^\x,t) + X^K^^\x,t) + ... (A.27) 
also converges uniformly then the solution of the integral equation (lA.ip . 



f{x)=u{s)-xj K{x,t)u{t)dt, (A.28) 
is represented by the 'reciprocal integral equation' 

u{x) = fix) + X I }C{x,t)f{t)dt. (A.29) 

J a 

The function /C(x, t) = }C{x, t; X) is therefore called the reciprocal or resolventkeinel. 

Series flA.26P and flA.27p are known as Neumann's series. They converge uni- 
formly if |A| is sufficiently small. Thus, for sufficiently small |A|, the reciprocal kernel 
is an analytic function of A. 

If the kernel K{x, t) is symmetric we can find a remarkably simple form for the 
resolvent kernel. Assuming again |A| to be sufficiently small, we make use of the 
series expansion (]A.20p for the symmetric kernels K^'^\x.,t).,K^^\x.,t)., ... and take 
the sum of the geometric series occurring in (1A.27P : we immediately obtain 

lC{x^ t- A) = K{x, t) + xf2 ^^l^M). (A.30) 

~l Ail^Ai - A) 

We see that the series on the right converges for every value of A which is not an 
eigenvalue, and the convergence is uniform in x and t. 

Relation (1A.30P which has been proved only for sufficiently small |A| provides 
the analytic continuation of the resolvent )C{x, t; A) into the entire complex A-plane, 
with the eigenvalues A, appearing as simple poles. Thus the resolvent of a symmetric 
kernel is a meromorphic function of X which possesses simple poles at the eigenvalues 
of the integral equation. Its residues at the poles Aj provide the eigenf unctions 
belonging to these eigenvalues. 



A. 3. Equations of the first kind 
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A. 3 Equations of the first kind 

Before trying to solve Fredholm integral equations of the first kind flA.4p we may first 
ask whether a solution does exist at all. We emphasise this point, since the theory 
is rather restrictive for the existence of solutions if compared to that of Fredholm 
equations of the second kind. 

The difficulty of finding a solution for the Fredholm integral equation (lA.4p 
stems from the fact that the integration operation over the sought solution u{t), 
is a smoothing process, especially, when combined with a nicely behaved kernel 
K{x,t). This means, for example, that if the solution u{x) is piecewise continuous, 
then the above integration operation on the r.h.s. of flA.4|) . with a continuous kernel 
K{x,t), would result in a smoother, i.e., continuous output f{x) of flA.4l) . as the 
given function at hand. So, if we are given a continuous function f{x), we cannot, 
in general guarantee an answer in the search for a solution u{x) among the class of 
continuous functions! In other words, if we look at the r.h.s. of (]A.4|) as an integral 



transform of piecewise continuous functions, then this transform maps such a class 
of functions to a more restrictive one, in this case continuous functions. Indeed, 
for more smooth kernels, i.e., if K{x,t) is differentiable, the class of piecewise or 
even integrable functions u{t) is mapped into differentiable functions f{x). Hence, 
for a continuous kernel K{x,t) (in both x and t) and a continuous output f{x), the 
integral equation of the first kind flA.4|) cannot, in general, be solved by a continuous 



function u{t). Of course, if K{x,t) is not very regular, then it is possible that this 
irregularity is combined with the smoothness of the integration operation and a 
continuous solution u{t), to produce a continuous output f{x). 

Theorem A. 3 For a continuous real and symmetric kernel, and a continuous f{x), 



a solution to (A. 4) exists only if the given function f{x) can be expressed in a series 



of the eig en functions of the kernel K{x^t), i.e., only if 

oo „fc 

f{x)=^ak(pk{x), ak= f{x)ipk{x)dx, (A.31) 
k=i -^^ 

where we have used the orthonormal set of eig enf unctions {ipk{x)} on [a,h]. □ 
With the condition stated in this theorem, the solution of flA.4p takes the form 

oo 

^\kak^k{x). (A.32) 



U{X) 

k=l 



We may note here that, although the existence of the (continuous) solution of 
flA.4|) in the form of the series flA.32|) is guaranteed, it is by no means a unique 



solution. This is the case, since if we add to the series in flA.32|) a function iIj{x) 
that is orthogonal to the kernel K{x,t), i.e., 

f K{x,t)ij{t)dt = 0, (A.33) 

J a 
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and substitute in flA.4|) . we obtain the same output f{x). So, for a unique solution 



u{x) in (lA.32p . we must insist that there are no functions iplx) that are orthogonal 
to the symmetric kernel K{x,t). 

The general case is covered by a theorem of Picard [PilO] which gives necessary 
and sufficient conditions for the existence of a square integrable solution u{t) of an 
equation of the first kind with an arbitrary (possibly unsymmetric) kernel. 

It is also useful to remark that Fredholm integral equations of the first kind can 
easily show the characteristics of ill-posedness defined in Section 11.31 We have just 
established a unique solution for Fredholm equations of the first kind for the specific 
case of symmetric kernels. What remains, for the present discussion, is the question 
whether the third quality, the stability of the solution of the problem, can also be 
established. Unfortunately, from the solution u{x) in (]A.32p we can show that the 
problem is not stable. When perturbing the given data /(x) by a small Sf{x) we 
expect the solution u{x), in its series representation (lA.32p . to be perturbed by a 
series representation of Sf{x) (or a constant multiple of it). However, one finds a 
magnification, i.e., a much larger corresponding change 6u in u{x). This is due to 
the eigenvalues factor in flA.321) . where they are increasing. If we write flA.321) . 
using the explicit form for a^, we have 



U[X} 

k=l 



Y^XkMx) / fiy)My)dy- (A.34) 
k=i -^^ 

Now, if we perturb /(x) by Sf{x) we get 

oo „5 

u{x) + 6u{x) = ^XkVkix) [f{y) + Sf{y)]fk{y)dy 



k=l 



^XkVkix) I fiy)(Pkiy)dy + ^Xk(pkix) j 6f{y)(pkiy)dy, 

(A.35) 



k=l k=l 



oo „b 

6u{x) = ^Xkekipk{x), ek= 6f{y)(pk{y)dy, (A.36) 
k=i -^^ 

where is the above expansion coefficient of the small perturbation 6f{x). In 
flA.361) we see that the expansion coefficients corresponding to the perturbation Su 
are magnified by a multiplicative factor A^, increasing with k, which clearly will 
cause 6u{x) to be a larger change compared to the given small change 6f{x). Hence, 
it seems that the solution u{x) is not stable in the Fredholm equation of the first 
kind (lA.4p . This possible ill-posedness adds to other difficulties for the existence of 
solutions for problems with more general kernels. It is no surprise that Fredholm 
integral equations of the first kind are in the forefront with regard to the research 
priority in the field of integral equations. 



Appendix B 

Green's function 



Green's functions provide a useful tool for dealing with source terms in differential 
equations. They are named after the British mathematician George Green, who 
first developed the concept in the 1830s. In this short overview we will present 
the general theory concerning the ordinary differential equations. The theory for 
partial differential equations is similar and it will not be treated here. Instead, some 
examples of the Green's function for the Laplace equation are presented. Finally, 
the Neumann Green's function for the unit disc is discussed in more detail. An 
elaborate theory of Green's functions can be found for example in |CH53a[ Wy28| . 

Throughout the present work, the method of Green's function was used to trans- 
form the partial differential equation governing the EIT into an integral equation 
(Chapters [9] and [T0|) . Furthermore, the Neumann Green's function is needed in 
Chapter [10] to solve these equations for the particular case of the unit disc. 



B.l General theory 



Let us consider a linear homogeneous self-adjoint differential expression of second 
order 



Lu = — 

ax 



du 
^ dx 



(B.l) 



for the function u{x) in the fundamental domain Q: a < x < b, where p, dp/dx 
and q are continuous functions of x and p > 0. The associated non- homogeneous 
differential equation is of the form 



Lu = f{x) 



(B.2) 



where (p{x) denotes a piecewise continuous function defined in Q. We are concerned 
with the boundary value problem: Find a solution of Eq. ( li?.^) which satisfies given 
homogeneous boundary conditions at the boundary of VL. 
For L we have the Green's theorem identity: 



.6 

/ dx [vLu — uLv] 

J a 



du 



dv\ 
dx J 



(B.3) 



In order to solve Eq. flB.2p for a general source (p{x) we introduce the so called Green's 
function G{x,x'), which is the solution for a point source at x': 



LG{x, x) = 6{x — x'), 



(B.4) 
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5{x — x') being the Dirac delta function. We then use the hnear character of (IB.ip 
to find the solution for a source v?(x) by superposition: 



dx'G{x, x'){p{x'). 



(B.5) 



In order to elaborate on the previous remarks, with suitable attention to the 
boundary conditions, let us suppose that the solution of (1B.4I) can be found. First, 
let us apply Green's theorem flB.3l) to u and v = G: 



dx [G{x, x')Lu{x) — u{x)LG{x, x')] 

du{x 



p{x) I G{x, x')- 



d 



— u{x)—G{x,x') 
dx dx 



-\ b 



(B.6) 



dxG{x, x'){p{x) — u{x) 



where we have used (IB. 21) . (IB. 41) and the integral properties of the delta function. 
Interchanging x and x' we find 



U[X] 



dx'G{x', x)ip{x')- 



pix) ( Gix'^x) —u{x)——;G{x\x)\ 
dx dx J 



■ (B.7) 



We can now consider several possibilities for the boundary conditions on u{x) at 
X = a,b. 

1. u{a) and u{b) given. 

For this case choose as boundary conditions on G{x',x) 

G{a,x) =G{b,x) = 0. (B.8) 



This has the desirable advantage of eliminating from Eq. (1B.7P the unknown quanti- 
ties du/dx' at x' = a, b. In fact (1B.7P reduces to 



u{x) 



dx'G{x', x)ip{x') + 



p{x')u{x') - — G{x' 
dx 



XI 



:b.9) 



This is the solution of Eq. (lB.2p in terms of known quantities. 
2. u{a) and u'{b) given. 

For this case choose as boundary conditions on G{x',x) 



d 



G(a,x) = 0, -—G(x', 
dx' 



X 



0. 



(B.IO) 



x'=b 



B.l. General theory 
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Again the unknown quantities, in this case u{b) and du/dx aX x = a are ehminated 
form Eq. (1B.7P and we obtain the solution of (]B.2p in terms of known quantities: 



U[X] 



dx'Gix' , x)(p(x') — p(a)u(a) - — G(x',x) 

dx' 



p{b)G{b,x)^ 



x'=b 

B.ll) 



3. Au{a) + Bu'{a) = X and Cu{b) + Du'{b) = Y given. 
Choose 



AG(a, x) + B —Gix, x) 
dx 



GG(b,x)+D^G(x',x) 
dx 



0, 



0. 



x'=b 



The quantity needed in flB.7p at x' = a is then 



G(a, x] 



du{x') 



dx' 



u(a) -^G(x',x) 
dx 



G{a, x) 

G{a, x) 
B 

G{a, x) 
B 

1 d 
A dx 



du{x') 



dx' 



— u(a] 



- — G{a,x) 



Au{a) + B 



du{x') 



dx' 



X 



-G{x' , x) 



X, 



(B.12) 



fB.13) 



which is a known quantity. A similar manipulation shows that the quantity needed 
at x' = 6 in (1B.7P can be expressed in terms of the given quantity Y . 

In each case we can choose the boundary conditions on G so as to eliminate 
unknown values of u and u' at the boundaries. In this way we obtain a solution 
of the inhomogeneous differential equation (1B.2|) once we find the Green's function 
G{x',x). 

We can prove that the Green's function is symmetric in its two arguments. To 
see this, apply the Green's theorem (IB.SP to u{x) = G{x,x') and v{x) = G{x,x"). 
Using Eq. (lB.4p we find 



G{x',x")-G{x",x') 



d d 
p(x) ( G(x, x")—G(x, x') — G(x, x')—G(x, x") 
dx dx 



. (B.14) 
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For the general case with boundary conditions (IB. 121) the r.h.s. of Eg. (IB. 141) vanishes 
and we find 



G{x',x") = Gix",x'). 
Let us now return to the differential equation flB.4p for G{x, x') 



d_ 

dx 



p(x)—G{x, x') 
dx 



q{x)G{x, x) = 6{x — x). 



(B.15) 



(B.16) 



For X ^ x' this reduces to the homogeneous equation 



d_ 

dx 



p{x)-^G{x, x') 
dx 



q{x)G{x,x') = 0, Xy^x', 



(B.17) 



so we can use whatever techniques are available for homogeneous equations in solving 
for G{x, x') in the two regions x > x' and x < x' . At a; = x' we must match the two 
solutions so obtained. The delta function 6{x — x') on the r.h.s. of flB.16|) is to be 
interpreted as due to a discontinuity in dG{x,x')/dx at x = x'. In fact, integrating 
flB.161) over an interval a; = x' — e to x = x' + e of infinitesimal length about x = x' 
and assuming q{x)G{x, x') is finite in this region so that 



x'+e 



dx q{x)G{x,x') — > 0, 

£— >0 



(B.18) 



we find 



^G{x,x'] 
dx 



p{x') 



:b.i9) 



provided p{x) is continuous at x = x'. On the other hand, we must take G{x, x') to 
be continuous at x = x': 



G{x' + e,x') = G{x' -e,x'). 



(B.20) 



If there were a discontinuity in G(x, x') itself, the first derivative dG{x, x')/dx would 
contain a term proportional to 5{x — x') and the second derivative d'^G{x, x')/dx'^ a 
term proportional to d6{x — x')/dx, in contradiction to the defining Eq.( ]B.16"l) . 

The homogeneous Eq. (lB.17p in the two regions x < x' and x > x' together with 
the two conditions (1B.19P and (1B.20P at x = x' and the boundary conditions at 
x = a and x = b discussed previously, completely determine the Green's function, 
the four conditions fixing the two pairs of integration constants obtained on solving 
the differential equation in the two regions. 

In the case of partial differential equations with homogeneous boundary con- 
ditions. Green's function can again be introduced as the kernel of an equivalent 
integral equation with the same properties as in the case of the ordinary differential 
equations |CH53a] . We will not repeat the theory, but rather give some examples. 
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B.2 Examples of Green's function for partial di- 
fferential equations 



Let us consider the Laplace equation for a two-dimensional domain fl 
A(f){x) = 0, with A = — + 



(B.21) 



the Laplace operator. In the following we will consider a few special cases: 

• Choose Go to be the fundamental (free) Green's function, satisfying: 

A^Go(x, x') = -5(x - x') for x, x' e Q. (B.22) 

This function is independent of the shape of Q and for a two-dimensional space 
it is: 



Go{x, x') = -— log \x - x'\ . 

ZTT 

Define Gd to be the Dirichlet Green's function given by 

A.xGd{x, x') — —5{x — x') for x,x' 
Gd{x, x') = for a; e and x' 



(B.23) 



(B.24) 



In contrast to Go, Gd depends on the shape of Q, and the difference between 
them is a harmonic term. For a disc of radius Gd is given by the formula 



+ r'^ - 2rr' cos{0 - 0') 
R'^ + ^r'^ - 2rr' cos{e - 6') 



(B.25) 



where (r, 9) and (r', 9') are the polar coordinates associated to x and x' re- 
spectively. 

Now consider Gat to be the Neumann Green's function which obeys the con- 
ditions: 



A^Gn{x,x') 

dn 



-5{x — x') + — for x,x' 



(B.26) 



{X, X 







for cc e and a;' e Q 



where \VL\ is the volume of the domain Vt. Like G/), Gjv also depends on the 
shape of the domain Vt and differs from Go by a harmonic term. For a disc of 
radius R, Gn reads: 



GN{r,9;r',9')^ -—log 



2 I '2 

r + r 



R' 



,2 



i?2 



2rr' cos{9 - 9') 
r'^ - 2rr' cos(9 - 9') 



(B.27) 



where by (r, 9) and (r', 9') we denoted again the polar coordinates associated 
to X and x' respectively. 
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B.3 Neumann Green's function for the unit disc 

The eigenvalue problem 



Av + Xv = 0, 
dv 



dr 



0, 



(B.28) 



r=l 



defines the eigenfunctions Vn of the Neumann Green's function, Eq. (]B.27p . assumed 
to have the eigenvalue expansion: 



Gn{x,x') = ^ — 



Vr,(X]V„(X 



The Laplacian in polar coordinates reads: 



(B.29) 



A 



d f d 



+ 



d fl d 



dr \ dr J 36 \r 36 
Let us search for solutions of flB.28p of the form f (r, 6) = f{r)u{6): 



(B.30) 



+ /^M = 



r dr dr"^ 



+ A- 



/ = 0, 



dl_ 

dr 



(B.31^ 



The solutions of the first equation are just the well-known eigenfunctions for the 
unit circle 



1 



sinZ^ for j = 1 
cos 16 for j = 2 



(B.32) 



'vr(l + 5^o) 

To solve the second equation in flB.311) . we make the change of variable: 

kr = p with A = k"^, (B.33) 
and find Bessel's equation 



P 



-/(p) + r(p)+ 1 



f{p) = 0, 



which has the regular solution 
f{p) = Mp) = Mkr), 



(B.34) 



(B.35) 



where Ji are the Bessel functions. Imposing the boundary condition (see Eq. (IB.26]) ). 
one finds 



kJ'iik) = 



(B.36) 



B.3. Neumann Green's function for the unit disc 
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and thus k are the zeros of the first derivative of Ji, usually denoted jf^. The 
eigenfunctions corresponding to the eigenvalues A = (j^m)^ will then be 

fiir) = JiULr). (B.37) 

One still has to normalise them. Recalling that 

I' JiULr)dr = ^(^l-l^j Jf{j*j, (B.38) 

the normalisation constants become 



Finally, the eigenfunctions corresponding to the eigenvalues A = (jf^)^, defined 
by the boundary eigenvalue problem flB.28p . are: 

or shortly 

viJr, e) = Q^Mj*^r)uiie), I e N+, (B.41) 
and the Neumann Green's function is 

G^(r, 6; r', = £ E "'-(^' ^-^^ ' (B.42) 

1,171=1 J = l 

On the boundary, i.e., on the unit circle, one has 

OO 2 CXD „ 

G^(l, e- 1, ^0 = E E E -^T^2 ■ (B-43) 

1=1 j=l m=l "''"^ 

Recalling that: 



m 



one finds: 

OO 2 

G^(l, 1, ^0 = E E (B.45) 

«=i j=i 



Appendix C 

SVD of an integral operator 



One of the most fruitful tools in the theory of linear inverse problems is the sin- 
gular value decomposition of a matrix and its extension to certain classes of linear 
operators. As we have seen in Chapter [21 SVD is basic both for understanding the 
ill-posedness of inverse problems and for describing the effect of the regularisation 
methods. Here, we will restrict the discussion to the SVD of integral operators with 
square-integrable kernel, needed in the EIT application of Section 110.11 
Let us define an integral operator of the following form: 

{Af){x)= I K{x,x')f{x')dx\ xen' (C.l) 
Jn 

where VL and Vt' are the object and the image domain respectively. 

We assume, for simplicity, that both the object and the image are square- 
integrable functions of the space variables, so that we have X = L^(fi) and 3^ = 
L^(f2'). Then Eq. flC.ip defines an operator from L'^{Q) into L'^{Q'). This operator 



is continuous if the integral kernel is square-integrable, i.e. 

\\K\\^= [ dx [ dx'\K{x, x')]"^ <oo. (C.2) 
Jn' Jn 

An integral operator satisfying this condition is usually called an integral operator 
of the Hubert- Schmidt class. 

In order to show that this operator is continuous we apply the Schwartz inequality 
to the r.h.s. of Eq. (lC.ip which, for fixed x, results in the scalar product of two square 
integrable functions. It follows 



\{Af){x)\'<(^ljKix,x')\'dx'^ (^jjf{x')\'dx' 



(C.3) 



If we integrate both sides of this inequality with respect to aj, and take the square 
root of the result, we get 

\\Af\\y<\\K\\.\\f\U, (C.4) 

i.e., the operator is bounded. This property implies the continuity of the operator. 
The adjoint A* of the operator A is given by 

{A*g){x')= [ K*{x,x')g{x)dx, x' e Q. (C.5) 
Jn' 
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C. SVD of an integral operator 



Then we can introduce the operators A = A* A and A = AA*. Both are integral 
operators with integral kernels given by 



A: K{x,x')= I K*{x",x)K{x",x')dx"; x,x'en, 



(C.6) 



A: k{x,x') = [ K{x,x")K*{x',x")dx"; x,x'En'. (C.7) 
Jn 

The integral operators A, A have the following properties: 

• Both operators are self-adjoint, i.e., for any pair of functions f,hmX and 
any pair of functions g, w in y 

(Af, h)x = if, Ah)^, (Ag, w)y = {g, Aw)y. (C.8) 
This property is a consequence of the following relations 

K*{x, x') = K{x', x), k*{x, x') = k{x', x) (C.9) 
which can be easily checked by means of Eqs. flC.6l) and flC.7l) : 

• Both operators are of the Hilbert-Schmidt class because their integral kernels 
are square-integrable (the proof of this result is similar to the proof of the 
inequality flC.4|) : the starting point is the application of the Schwartz inequality 
to the r.h.s. of Eqs.flOel) and (^Cl\l ): 



• Both operators are positive semi-definite 

{Af,f)x>0, {Ag,g)y>0. (C.IO) 

Remark C.l We sketch the proof of this property in the case of A. Indeed, from 
Eq. liC. hy means of an exchange of the integration order we have 



{Af, f) 



X 



[ ( [ Kix,x')fix')dx') f*{x)dx 



K{x" , x)f(x)dx 



dx" > 0. (C.ll) 



A similar proof applies to the case of A. □ 

According to the Hilbert-Schmidt theory [Mi57j . an integral operator with a 
symmetric and square-integrable kernel has real eigenvalues with finite multiplicity. 
Moreover the eigenfunctions associated with different eigenvalues are orthogonal. 
The eigenvalues form, in general, a countable set with zero as accumulation point. 

If IC{x,x') is a symmetric and square-integrable kernel, let Ai,A2,A3,... be the 
sequence of the eigenvalues of the corresponding integral operator, ordered in such 
a way that |Ai| > IA2I > IA3I > each eigenvalue being counted as many times 
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as its multiplicity. Moreover, let Vi{x), V2(x),V3(x), ... be the sequence of the eigen- 
functions associated with these eigenvalues. They constitute an orthonormal set of 
square-integrable functions. Then the basic result of the Hilbert-Schmidt theory is 
the following spectral representation of the kernel )C{x, x') 

oo 

}C{x,x') = J^XkV,{x)vl{x'), (C.12) 

k=l 

the series being convergent in the sense of the L^-norm. Accordingly, the eigenvalues 
Afc satisfy the following condition 

oo 

^A^<oo. (C.13) 

k=l 

The result stated above apply to the operators A and A, whose non-zero ei- 
genvalues are positive because they are positive semi-definite operators. One can 
also show that the operators A and A have the same eigenvalues with the same 
multiplicity. 

Let us denote by al the positive eigenvalues of A, with the ordering o"^ > (t| > 
■■■ > CTfe > ••• and cr^ for A; cxD. Again each eigenvalue has been counted as 
many times as its multiplicity, which is certainly finite as follows from the Hilbert- 
Schmidt theory. 

A normalised cigenfunction is associated to each eigenvalue o"^ and these 
eigenfunctions constitute an orthonormal system, i.e. {vk,Vj)x = ^kj- Then, to each 
cigenfunction Vk of A we can associate a function Uk in y defined by 

Uk{x)^—{Avk){,x). (C.14) 

Using the relation AA — AA, which can be easily proved by an exchange of integra- 
tion order in the definition of the integral kernels, we obtain 



{Auk){x) = —{AAvk){x) = ai —{Avk){x) = aiuk{x) (C.15) 

\(^k 



and also 



{Uk,Uj)y = -^{Avk,Av.j)y 

= -^{vk,Avj)x ^—{vk,Vj)x = hj- (C.16) 

Therefore all eigenvalues a1 of A are also eigenvalues of A and the are the co- 
rresponding eigenfunctions, which constitute an orthonormal system in y . In order 
to show that in this way we have obtained all eigenvalues and eigenfunctions of A 
it is sufficient to repeat the same argument starting from A. If its eigenvalues and 
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C. SVD of an integral operator 



eigenfunctions are denoted by ak and Uk respectively, then we can show, by means 
of the relation AA* = A* A, that the functions of X defined by 

Vk{x) = —{A*uk){x) (C.17) 

Cfk 

are orthonormal eigenfunctions of A associated with the eigenvalues a^. 

Eqs. (IC.14|) and flC.171) define the usual shifted eigenvalue problem, which we 
write now explicitly in terms of the integral operators 

/ K{x,x')vk{x')dx' = akUk{x), (C.18) 
Jq 

/ K* {x,x')uk{x)dx = akVk{x'). (C.19) 

The proof of the SVD of the operators A and A* requires the completeness of 
the Hilbert spaces X and y. We do not report this proof here (see, for instance, 
|Gr93] ). We only give the result, which takes the form: 

oo 

{Af){x) = Y,^kU.Vk)xUk{x) (C.20) 
fc=i 

and also 

oo 

{A*f){x) = Y,^k{9^Uk)yVk{x). (C.21) 

k=l 

the series being convergent with respect to the norms of y and X, respectively. 

We also remark that the SVD of A is equivalent to the following series expansion 
of the kernel K{x, x') 

oo 

K{x, x') = (rkUk{x)vl{x'), (C.22) 

k=l 

the series being convergent with respect to the L^-norm. From this expansion (which 
is a generalisation of equation (1C.12P ). from the orthogonality of the singular func- 
tions and from the definition (lC.2p of H-ft'lP it follows that 

oo 
k=l 

and therefore the sum of squares of the singular values of a Hilbert- Schmidt integral 
operator is convergent. 

Using the SVD (]C.20p it is also possible to show that the norm of the integral 
operator A is given by 

\\A\\ = a^. (C.24) 

From this equation and (lC.23p . it follows that the norm of the integral operator is 
never greater than the L'^-norm of the integral kernel; the two norms coincide if and 
only if the integral operator has rank one. 
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Statistics 



In statistics we are interested in using a given sample of data to make inferences 
about a probabilistic model, e.g., to assess the model's validity or to determine the 
values of its parameters. There are two main approaches to statistical inference, 
which we may call frequentist and Bayesian. In frequentist statistics, probability is 
interpreted as the frequency of the outcome of a repeatable experiment. Note that 
in frequentist statistics one does not define a probability for a hypothesis or for a 
parameter. 

Frequentist statistics provides the usual tools for reporting objectively the out- 
come of an experiment without needing to incorporate prior beliefs concerning the 
parameter being measured or the theory being tested. 

In Bayesian statistics however, the interpretation of probability is more general 
and includes the notion of a "degree of belief". One can then speak of a proba- 
bility distribution function (p.d.f.) for a parameter, which expresses one's state of 
knowledge about where its true value lies. Bayesian methods allow for a natural 
way to include additional information such as physical boundaries and subjective 
information; in fact they require as input the ^^nor p.d.f. for the parameters, i.e., the 
degree of belief about the parameters' values before carrying out the measurement. 

Bayesian techniques are often used to treat systematic uncertainties, where the 
subjective beliefs about, say, the accuracy of the measuring device may enter. 
Bayesian statistics also provides a useful framework for discussing the validity of 
different theoretical interpretations of the data. 

For many inference problems, the frequentist and Bayesian approaches give the 
same numerical answers, even though they are based on fundamentally different 
interpretations of probability. For small data samples, however, and for measure- 
ments of a parameter near a physical boundary, the different approaches may yield 
different results, so we are forced to make a choice. For a discussion of Bayesian vs. 
non-Bayesian methods, see references written by a statistician |Ef86] . by a physicist 
[ToQS] . or the more detailed comparison in Ref. |50A98] . 

In what follows we will concentrate on frequentist statistics and briefly discuss 
least-squares parameter estimation, confidence regions and the t^st of goodness- 
of-fit. In Chapters [5l [6] and [7] we aimed to determine values of some fundamental 
QCD parameters, the so called condensates, in a least-squares sense. The errors on 
these parameters were then defined through constant boundaries as CLs around 
the values found. It was also important to define the level of agreement between 
data and theory, for which we have made use of the test. 
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D. Statistics 



D.l Least-squares parameter estimation 

Consider a situation in which N data pairs, Xi, tji, with i = 1,2,.. .,N, are to be 
modeled by a functional relationship 

F{x)=F{x,e), y, = F{x„e), (D.l) 

where 6 is the vector of M adjustable parameters 6j. The goal is to determine the 
"best" values of 6 that describe the true values ^o- 

Given a particular data set of Xi's and yi's, we have the intuitive feeling that 
some parameters are very unlikely (those for which the model function F{x) looks 
nothing like the data) while others may be very likely (those that closely resemble 
the data). How can we quantify this intuitive feeling? How can we select fitted 
parameters that are most likely to be correct? It is not meaningful to ask the 
question, "What is the probability that a particular set of fitted parameters is 
correct?" The reason is that there is no statistical universe of models from which 
the parameters are drawn. There is just one model, the correct one, and a statistical 
universe of data sets that are drawn from it. 

That being the case, we can, however, turn the question around, and ask, "Given 
a particular set of parameters, what is the probability that this data set could have 
occurred?" If the ?/j's take on continuous values, the probability will always be zero 
unless we add the phrase, "...plus or minus some fixed Ay on each data point". 
So let's always take this phrase as understood. If the probability of obtaining the 
data set is infinitesimally small, then we can conclude that the parameters under 
consideration are unlikely to be right. Conversely, our intuition tells us that the 
data set should not be too improbable for the correct choice of parameters. 

In other words, we identify the probability of the data given the parameters as 
the likelihood of the parameters given the data. This identification is entirely based 
on intuition. It has no formal mathematical basis in and of itself. Once we made 
this intuitive identification, however, it is only a small further step to decide to fit 
for the parameters precisely by finding those values that maximise the likelihood 
defined in the above way. 

Suppose that each data point yi has a measurement error ai and that yi are 
independently random and distributed as a normal (Gaussian) distribution around 
the true model F{x). Then the probability P that the data set is described by F{x) 
is given by 

-luPocx'W = g( ''--^|"--^> ) . (D,2) 

Maximising the joint probability function, P, is equivalent to minimising x^- 
This is achieved by differentiating Eq. flD.2p for each 6j 



V _ f m-F{x,,0) \ ( dF{x,,0) \ 



D.2. Confidence regions 



171 



and equating to zero. Obtaining Eq. (]D.3P for eacli j yields a set of M equations, 
one for each unknown 6j. Minimising thus reduces to solving a system of M 
equations. 

The curvature matrix 

N 



lav 



2 86.6 



1=1 ' 



dF{xf, 6) dF{xf, 0) 



36, 



06, 



defines the covariance matrix for the parameters: 



V 



a 



-1 



(D.4) 



(D.5) 



Once is minimised, the variances in the parameters 6 are approximated by the 
diagonal matrix elements of V: 

J2 



(7a 



(D.6) 



D.2 Confidence regions 

The full probability distribution is a function defined on the M-dimensional space 
of parameters 6. A confidence region (or confidence interval) is just a region of that 
M-dimensional space that contains a certain percentage of the total probability 
distribution. When pointing to a confidence region one can say, e.g., "this is a 99% 
chance that the measured values Ui are obtained from the model with the parameter 
values 6 lying within this region" . 

It is worth emphasising that one has to determine both the confidence level 
(99% in the above example), and the shape of the confidence region. The only 
requirement is that the region does include the stated percentage of probability. 
Certain percentages are, however, customary in scientific usage: 68.27% (the lowest 
confidence worthy of quoting) also called laCL, 90%, 95.45% (2crCL), 99% and 
99.73% (ScxCL). Higher CL's are conventionally "ninety-nine point nine ... nine". 
As for shape, obviously one wants a region that is compact and reasonably centred 
on since the whole purpose of a confidence limit is to inspire confidence in 
that value. In one dimension, the convention is to use a line segment centred on 
the measured value; in higher dimensions, ellipses or ellipsoids are most frequently 
used. 

The numbers 68.27%, 95.45% and 99.73% and the use of ellipsoids are connected 
with a Gaussian distribution, however not allways relevant. In general, the proba- 
bility distribution of the parameters will not be Gaussian, and the above numbers, 
used as levels of confidence, are purely matters of convention. In these cases, the 
simple quotation of error ranges like i do not provide information about the 
shape of the p.d.f. 

Fig jD.ll sketches a possible probability distribution for the case M = 2. Shown 
are three different confidence regions, all at the same confidence level. The two 
vertical lines enclose a band (horizontal interval) which represents the luCR for the 
variable 6i without regard to the value of 62- Similarly the horizontal lines enclose 
a IcrCR for 62- The ellipse shows a joint laCR for pairs of values 61 and 62. 
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Figure D.l: Confidence intervals in 1 and 2 dimesions. The same fraction 
of measured points (here 68.27%) lies (i) between the two vertical lines, (ii) 
between the two horizontal lines, (iii) within the ellipse. 



D.3 Constant boundaries as CL 

When the method used to estimate the parameters is x^-minimisation, as in Section 
ID. 11 then there is a natural choice for the shape of confidence intervals, whose use 
is almost universal. For the observed data, the value of at is a minimum. Call 
this minimum Xmin- If the vector 6 of parameter values is perturbed away from Oq, 
then increases. The region within which increases by no more than a chosen 
amount A^^ defines some M-dimensional confidence region around ^o- If is 
set to be a large number, this will be a big region; if it is small, it will be small. 
Somewhere in between there will be choices of Ax^ that cause the region to contain, 
variously, la, 2cr, etc. of probability distributions for 0, as defined above. These 
regions are taken as the confidence regions for the parameters ^o- Values of Ax^ for 
different confidence regions and different numbers of parameters are listed in Table 

EH 

Very frequently one is interested not in the full M-dimensional confidence region, 
but in individual confidence regions for some smaller number v of parameters. For 
example, one might be intersted in the confidence interval of each parameter taken 
separately (the bands in Fig jD.l|) . in which case u = 1. In that case, the natural 



confidence regions in the //-dimensional subspace of the M-dimensional parameter 
space are the projections of the M-dimensional regions defined by fixed Ax^ into 
the i/-dimensional spaces of interest. 
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CR(%) 


M = 1 


M = 2 


M = 3 


68.27 


1.00 


2.30 


3.53 


90.00 


2.71 


4.61 


6.25 


95.45 


4.00 


6.18 


8.03 


99.00 


6.63 


9.21 


11.34 


99.73 


9.00 


11.83 


14.16 



Table D.l: Ax^ corresponding to a confidence region CR, for joint estimation 
of M parameters based on Gaussian (normal) p.d.f. 

D.4 The test 

'Don't ask what it means, but ruther how it is used. ' 

L. Wittgenstein 



Often one wants to quantify the level of agreement between the data and a 
hypothesis without explicit reference to alternative hypotheses. This can, in many 
cases, be done by the test invented in 1900 by Karl Pearson [PelQGG] . 

The goodness-of-fit is quantified by giving the observed significance level also 
called p-value. The p-value is defined as the probability to find in the region 
of equal or smaller compatibility with the considered hypothesis than the level of 
compatibility observed with the actual data. Since is defined such that large 
values correspond to poor agreement with the hypothesis, the p-value will be 

/■oo 

p= / f{z;,n)dz (D.7) 

where is the value obtained in the minimisation process (e.g., in the present 
work the value of ximin' Eq- (I5.10I) . was used to estimate the agreement between 
experimental data and the theoretical model, i.e., QCD). f{z]n) is the x^ p.d.f. 
with n degrees of freedom^, 

f(z;n) = _-l__zW2-ig-./2 P 8) 

■'^ ' ' r(n/2)2"/2 ' ^ ' 

where ^ > and f{z; n) = for z < 0. Here F denotes the Gamma function. 

Values of p can be obtained from Fig JD.2[ If the conditions for using the x^ 
p.d.f. do not hold, the p-value can still be defined as before, but the p.d.f., /, must 
be determined by other means, e.g., using a Monte Carlo calculation. 



^The term degrees of freedom is a measure of the number of independent pieces of information on 
which the precision of a parameter estimate is based. The degrees of freedom for an estimate equals 
the number of observations (values) minus the number of additional parameters estimated for that 
calculation. As we have to estimate more parameters, the degrees of freedom available decreases. It 
can also be thought of as the number of observations (values) which are freely available to vary given 
the additional parameters estimated. 
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If one finds a value much greater than n and a correspondingly small p- 
value, one may be tempted to expect a high degree of uncertainty for any fitted 
parameter. Although this may be true for systematic errors in the parameters, it is 
not in general the case for statistical uncertainties. If, for example, the error bars 
(or covariance matrix) used in constructing the ^"^^ underestimated, then this 
will lead to underestimated statistical errors for the fitted parameters. But in such 
a case an estimate 6q can differ from the true value by an amount much greater 
than its estimated statistical error. The standard deviations of estimators that one 
finds reflect how widely the estimates would be distributed as if one were to repeat 
the measurement many times, assuming that the measurement errors used in are 
also correct. They do not include the systematic error which may result from an 
incorrect hypothesis or incorrectly estimated measurements errors in the x^- 

Since the mean of the x^ distribution is equal to n, one expects in a "reasonable" 
experiment to obtain x^ k, n. Hence the quantity x^ /n is sometimes reported. 
Since the probability distribution function of x^ /n depends on n, however, one must 
report n as well in order to make a meaningful statement. The p- values obtained 
for different values of x^/n are shown in Fig JD.3[ 




Figure D.2: The p-value as a function of x^ foi" ^ degrees of freedom. 
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Degrees of freedom n 

Figure D.3: The "reduced" x^, equal to jn, for n degrees of freedom. 
The curves show jn that corresponds to a given p-value as a function of n. 
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